# Methods and Conversations in (Post)Modern Thermodynamics

Francesco Avanzini, Massimo Bilancioni, Vasco Cavina, Sara Dal Cengio, Massimiliano Esposito, Gianmaria Falasco, Danilo Forastiere, Jose Nahuel Freitas, Alberto Garilli, Pedro Harunari, et al.

► To cite this version:

Francesco Avanzini, Massimo Bilancioni, Vasco Cavina, Sara Dal Cengio, Massimiliano Esposito, et al. Methods and Conversations in (Post)Modern Thermodynamics. SciPost Physics Lecture Notes, 2024, pp.80. 10.21468/SciPostPhysLectNotes.80. [hal-04308373](https://hal.science/hal-04308373v1)

HAL Id: hal-04308373

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

Submitted on 27 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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 CC BY 4.0 - Attribution - International License
Methods and Conversations in (Post)Modern
Thermodynamics

Francesco Avanzini<sup>1,2,§</sup>, Massimo Bilancioni<sup>1,¶</sup>, Vasco Cavina<sup>1,©</sup>, Sara Dal Cengio<sup>3,•</sup>, Massimiliano Esposito<sup>1,†</sup>, Gianmaria Falasco<sup>4,ª</sup>, Danilo Forastiere<sup>4,5,§</sup>, Nahuel Freitas<sup>6,£</sup>, Alberto Garilli<sup>1,‡</sup>, Pedro E. Harunari<sup>1,§</sup>, Vivien Lecomte<sup>3,*</sup>, Alexandre Lazarescu<sup>7,∫</sup>, Shesha G. Marehalli Srinivas<sup>1,Nº</sup>, Charles Moslonka<sup>8,™</sup>, Izaak Neri<sup>9,||</sup>, Emanuele Penocchio<sup>10,¤</sup>, William D. Piñeros<sup>1,∞</sup>, Matteo Polettini<sup>1,•</sup>, Adarsh Raghu<sup>9,*</sup>, Paul Raux<sup>11,12,€</sup>, Ken Sekimoto<sup>8,13,Ξ</sup>, and Ariane Soret<sup>1,√</sup>

1 Department of Physics and Materials Science, University of Luxembourg, Campus Limpertsberg, 162a avenue de la Faïencerie, L-1511 Luxembourg (G. D. Luxembourg)

2 Department of Chemical Sciences, University of Padova, Via F. Marzolo, 1, I-35131 Padova, Italy

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

4 Department of Physics and Astronomy, University of Padova, Via Marzolo 8, I-35131 Padova, Italy

5 INFN, Sezione di Padova, via Marzolo 8, I-35131 Padova, Italy

6 Universidad de Buenos Aires, Facultad de Ciencias Exactas y Naturales, Departamento de Física. Buenos Aires, Argentina

7 Institut de Recherche en Mathématique et Physique, UCLouvain, Belgium

8 Laboratoire Gulliver, UMR CNRS 7083, ESPCI Paris, Université PSL

9 Department of Mathematics, King's College London, Strand, London, WC2R 2LS, UK

10 Department of Chemistry, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA

11 Université Paris Cité, CNRS, UMR 8236-LIED, 75013 Paris, France

12 Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France

13 Laboratoire Matière et Systèmes Complexes, UMR CNRS 7057, Université Paris Cité

$francesco.avanzini@unipd.it

¶massimo.bilancioni@uni.lu

©vasco.cavina@uni.lu

•sara.dal-cengio@univ-grenoble-alpes.fr

†massimiliano.esposito@uni.lu

ªgianmaria.falasco@unipd.it

§danilo.forastiere@unipd.it

£nfreitas@df.uba.ar

‡alberto.garilli@uni.lu

§§pedro.harunari@uni.lu

*vivien.lecomte@univ-grenoble-alpes.fr

⁷alexandre.lazarescu@gmail.com

№shesha.marehalli@uni.lu

™charles.moslonka@espci.psl.eu

||izaak.neri@kcl.ac.uk

¤emanuele penocchio@northwestern.edu

∞william.pineros@uni.lu

•matteo.polettini@uni.lu

*raghu.adarsh@kcl.ac.uk

€paul.raux@etu.u-paris.fr

Ξken.sekimoto@espci.psl.eu

√ariane.soret@uni.lu

November 3, 2023

Lecture notes after the doctoral school (Post)Modern Thermodynamics held at the University of Luxembourg, December 2022, 5-7, covering and advancing continuous-time Markov chains, network theory, stochastic thermodynamics, large deviations, deterministic and stochastic chemical reaction networks, metastability, martingales, quantum thermodynamics, and foundational issues.

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

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

Foreword 5

Introduction 5

Outline 7

1 Markov Chains and First-Passages 8

1.1 Introduction 8

1.2 Continuous-time Markov Chains: basic notions 8

1.3 Simulating a trajectory: Gillespie's algorithm 10

1.4 First-passage time problems 12

1.5 Other first event problems 17

1.6 "Markovianization" of non-Markov problems 19

1.7 Discrete-time Markov Chains 21

1.8 First-transition times and hidden networks 23

2 Network Thermodynamics 28

2.1 Introduction 28

2.2 Elements of graph theory 28

2.3 Network theory applied to CTMC: the Markov chain tree theorem 34

2.4 Network theory applied to observable decomposition 38

2.5 An example of chemical reaction network 40

3 Stochastic Thermodynamics 44

3.1 Stochastic dynamics 44

3.2 Stochastic thermodynamics for a single reservoir of energy 46

3.3 ST for multiple reservoirs 53

3.4 Internal entropy and matter transfer 54

3.5 Fluctuation theorems and thermodynamic uncertainty relations 56

4 Deterministic Chemical Reaction Networks 58

4.1 Introduction 58

4.2 Basic Setup 59

4.3 Dynamics 60

4.4 Local Detailed Balance 63

4.5 Thermodynamics of Closed CRNs 67

4.6 Thermodynamics of Open CRNs 70

4.7 Conclusions 78

5 Stochastic Chemical Reaction Networks 79

5.1 Introduction 79

5.2 Fishing for Poisson 79

5.3 Fishing for Poisson-like 83

5.4 Moieties and product-form 84

5.5 One trivial cycle 86

5.6 One less trivial cycle: Schlögl model 89

5.7 Zero deficiency 91

5.8 Deficiency and response 93

5.9 Considerations 98
<table><tbody><tr><td><b>6 Dynamical Large Deviations: at long times</b></td><td><b>99</b></td></tr><tr><td>6.1 Introduction: Why large deviations?</td><td>99</td></tr><tr><td>6.2 Basics and formalism</td><td>100</td></tr><tr><td>6.3 Application: Large Deviations in Diffusive Systems</td><td>102</td></tr><tr><td>6.4 Biased dynamics: making rare dynamics typical</td><td>105</td></tr><tr><td><b>7 Dynamical Large Deviations: at large size</b></td><td><b>110</b></td></tr><tr><td>7.1 Dynamical large deviations</td><td>110</td></tr><tr><td>7.2 Path integral formalisms in physics: a comparison</td><td>124</td></tr><tr><td><b>8 Metastability In and Out of Equilibrium</b></td><td><b>126</b></td></tr><tr><td>8.1 Introduction</td><td>126</td></tr><tr><td>8.2 Metastability for Markov processes</td><td>127</td></tr><tr><td>8.3 Overdamped Langevin dynamics: low temperature limit</td><td>129</td></tr><tr><td>8.4 Jump processes in occupation-number space: macroscopic limit</td><td>134</td></tr><tr><td>8.5 Markov jump process on the space of attractors</td><td>136</td></tr><tr><td>8.6 Conclusions</td><td>137</td></tr><tr><td><b>9 Martingale Approach for First-Passage Problems</b></td><td><b>138</b></td></tr><tr><td>9.1 Introduction</td><td>138</td></tr><tr><td>9.2 Gambler's ruin problem</td><td>139</td></tr><tr><td>9.3 Martingales: Definition, Represenations, and Properties</td><td>143</td></tr><tr><td>9.4 Martingale solution to the gambler's ruin problem</td><td>146</td></tr><tr><td>9.5 Thermodynamic versions of the gambler's ruin problem</td><td>148</td></tr><tr><td>9.6 Discussion</td><td>155</td></tr><tr><td>9.7 Appendix: Martingality of Eq. (9.81)</td><td>156</td></tr><tr><td><b>10 Quantum Thermodynamics</b></td><td><b>158</b></td></tr><tr><td>10.1 Introduction</td><td>158</td></tr><tr><td>10.2 Mathematical framework for quantum statistics</td><td>159</td></tr><tr><td>10.3 Quantum measurement and thermodynamic quantities</td><td>163</td></tr><tr><td>10.4 Quantum fluctuation theorems</td><td>167</td></tr><tr><td>10.5 Conclusion and perspectives</td><td>172</td></tr><tr><td><b>11 Why "(Post)Modern Thermodynamics"?</b></td><td><b>174</b></td></tr><tr><td>11.1 Introduction</td><td>174</td></tr><tr><td>11.2 Entry survey & breaking the ice</td><td>174</td></tr><tr><td>11.3 Discussion</td><td>176</td></tr><tr><td>11.4 Denoument</td><td>177</td></tr><tr><td>11.5 Considerations</td><td>179</td></tr><tr><td><b>A Proposed unified notation</b></td><td><b>181</b></td></tr><tr><td>A.1 Free symbols</td><td>181</td></tr><tr><td>A.2 Acronyms</td><td>181</td></tr><tr><td>A.3 Relations</td><td>182</td></tr><tr><td>A.4 Conventions</td><td>182</td></tr><tr><td>A.5 General</td><td>183</td></tr><tr><td>A.6 Linear algebra</td><td>183</td></tr><tr><td>A.7 Oriented graphs</td><td>184</td></tr><tr><td>A.8 Stochastic variables</td><td>184</td></tr><tr><td>A.9 Stochastic processes</td><td>184</td></tr></tbody></table>
<table><thead><tr><th>SciPost Physics Lecture Notes</th><th>Submission</th></tr></thead><tbody><tr><td>A.10 Markov chains</td><td>185</td></tr><tr><td>A.11 Thermodynamics</td><td>185</td></tr><tr><td>A.12 Chemical reaction networks</td><td>186</td></tr><tr><td><strong>References</strong></td><td><strong>186</strong></td></tr></tbody></table>
# Foreword

## Massimiliano Esposito

Progress in nonequilibrium physics has been quite spectacular over the last 25 years. The trigger was without doubt the discovery of fluctuation theorems in many seemingly different contexts (thermostatted dynamics, Hamiltonian dynamics, stochastic dynamics) and the need to come up with a coherent understanding of the relation amongst these various fluctuation theorems. This effort converged into what is often called nowadays, broadly speaking, stochastic thermodynamics, encompassing the description of both classical and quantum systems. But at the same time that stochastic thermodynamics was being consolidated, research in the field also diversified and started interfacing with many new research areas.

The doctoral school (Post)Modern Thermodynamics that was held at the University of Luxembourg in December 2022, 5-7, represents that evolution. The first three courses covered classical topics in stochastic thermodynamics: defining Markov jump processes, analyzing them with elements of network theory, constructing the thermodynamic quantities (e.g. entropy production, local detailed balance, affinities), and analysing the role of coarse-graining. The emphasis on first-passage times and martingales in courses 1 and 9 is already a move toward more recent developments in stochastic thermodynamics, which showed that various first-passage problems can be constrained by thermodynamics. Courses 6 and 4 reflect the importance that the study of chemical reaction networks using stochastic thermodynamics has taken over the last decade. These systems are ideal to study the fate of stochastic thermodynamics in the macroscopic limit. Indeed their deterministic description can be viewed as emerging from a large volume limit. Course 5 introduces large deviation theory, which provides a powerful tool to characterize how nonequilibrium fluctuations scale in such a macroscopic limit. Course 8 analyses instead metastability, another phenomenon arising in the macroscopic limit (or equivalently in the low noise limits), and uses it to provide key insight into the central concept of local detailed balance. Course 10 focuses instead on quantum thermodynamics, another field that has attracted significant attention over the last decade and that aims at extending the concepts from classical stochastic thermodynamics to quantum systems and making contact with concepts from quantum information theory.

Matteo Polettini took the lead in organizing this doctoral school and did it in a truly collective manner, involving many young researchers. He also added a philosophical flavor to it, with a session dedicated to discussing foundational issues in thermodynamics and maybe decoding the enigmatic title of the school. A report on this session is the object of the last course 11.

We received numerous very enthusiastic feedback from the participants showing that the event was a real success. The hope is that the present lecture notes will help many young researchers enter the rapidly evolving field of nonequilibrium thermodynamics.

# Introduction

## Pedro Harunari, Vasco Cavina, William Piñeros and Matteo Polettini

In early December 2022, amidst the first snowflakes and city lights illuminating Christmas markets, we gathered in the small city of Luxembourg for an immersive week dedicated to thermodynamics. The school+workshop provocatively named (post)Modern Thermodynamics brought together approximately 130 physicists to delve into ongoing research, established principles and interpretations of thermodynamics, and their recent unfoldings.
The first part of the event was a school aimed at, but not limited to, graduate students. Experts in their fields gave 2-hour-long lectures on topics such as stochastic thermodynamics, quantum thermodynamics, mathematical methods for statistical physics, and chemical reaction networks. Our goal was to provide introductory courses on relevant areas of research and improve the toolboxes of young researchers. More than textbook-style lectures, we asked lecturers to include a twist of their own research, making the process a bit more personal and communicating state-of-the-art developments to the students. This resulted in ten different lectures about the frontiers of thermodynamics by teachers with different original takes on a theory that puzzled physicists for a good part of the last three centuries.

The disparity of background, experiences and points of view of our teachers profoundly enriched the scientific discussion during the school but also compelled us to look for a way to increase its cohesion. In a manner befitting good physicists, we decided to run an experiment and to assign to each of our lecturers an “angel”, whose role was to help in preparing the lectures and avoid overlaps within the program. Another important task assigned to the angels was to help in standardizing the notation among the various contributions: Much confusion arises from assigning distinct names/symbols to the same ideas such as e.g. entropy production versus entropy production rate, order of indices, negative versus positive signs conventions for directed heat flows, etc. During the lectures and in this version of them, we and the angels did our best to dissipate this confusion and strived to present a treatment that is notationally and conceptually homogeneous.

Given the pedagogic character of the lectures, we prepared the following notes with the help of some young participants of the (post)Modern Thermodynamics school. The contribution of the students to the event was invaluable, their curiosity and excitement being the main driving forces behind all the scientific discussions. Having them collaborate on these lecture notes was therefore a natural way to try and capture their energy and enthusiasm on the exposition of these subjects.

Of course, a team composed primarily of students preparing the notes was not without its challenges despite the healthy number of volunteers. We identified room for improvement in distributing workload and coordinating the various perspectives from members of different seniorities. Nonetheless we found we could partially integrate these differences via an internal “peer-review” process where team members could discuss and receive feedback on each other’s content. Altogether, we believe the effort was an instructive exercise in collaborative academic writing for both students and organizers.

Our final communal experiment, and indeed the content of the last and namesake lecture 11 of the school, was an open-ended discussion to reflect on fundamental issues in and about thermodynamics. This took the form of a set of suggested questions, both by organizers and participants alike, in which everyone could also contribute an answer anonymously online for later in-person discussions. One of the most surprising outcomes of the survey was the lack of a common interpretation on the foundations of thermodynamics, showing how even in the presence of universally accepted experimental results and a well-established formalism there is still room for disparate conceptual beliefs. These differences sparkled an interesting philosophical debate that was highly appreciated by participants of the school+workshop and lead us to believe that there could be more space for this kind of discussion at scientific conferences.

The event itself was held free of charge following funds from already awarded grants. We were additionally able to offer financial aid, by means of free accommodation, to all students who requested it. This allowed participants from less fortunate backgrounds to join the school and attend in person. In total, 20 students benefitted from this opportunity, increasing the overall pool of students.

Post-(post)Modern Thermodynamics, we ran a survey that revealed that most partic-
ipants were satisfied with the chosen themes and overall organization. The criticism that
stood out was the lack of time for relaxation and discussion in view of the dense program.
We believe this problem could be circumvented by relaxing the number of lectures/talks
or by making clear that participants understand they are not obliged to attend all the
sessions as science is social, and conferences are made for exchanging.

We tried to make this school an opportunity for young researchers to grow, a “rich and challenging environment for the individual to explore”, using the word of Noam Chomsky (to whom our social event, a widely participated billiard tournament, was dedicated).

We hope these lecture notes and the reported experience will serve the community of thermodynamicists for many more such future events.

## Outline

The first three lectures, Chaps. 1, 2 and 3, deal with the dynamics of Markov processes defined over networks, and how thermodynamics is tied to them. Lectures 4 and 5 address chemical reaction networks from deterministic and stochastic standpoints, presenting both analytical and numerical approaches, with particular attention to their nonequilibrium thermodynamics. Lectures 6 and 7 address the large deviation principle for describing fluctuations of observables in the limits of long times and large systems, and Lecture 8 discusses in more detail the metastable states that arise between these limits. Lecture 9 introduces martingale theory and its recently discovered connections to thermodynamics. Lecture 10 overviews the mathematical framework of quantum statistics and presents a consistent approach to defining thermodynamic quantities at the quantum level. Finally, lecture 11 wraps up with a discussion on thermodynamics fundamentals considering the input of school attendees, as we elaborated above.

We prepared a tentative unified notation (see Appendix A) and invited lecturers/angels to discuss and employ it at will. By no means did we intend to impose a standard on the field, but rather implement a time-saving tool when reading across the different sets of lectures and their related concepts. Some lecturers agreed to it, some others preferred maintaining diversity as a pedagogical tool. Interestingly, the parts that were more agreed upon were the more standard mathematical ones (linear algebra, graph theory) while little agreement was found in the use of symbols for thermodynamic quantities, apart from S for entropy. In the meantime, the issue became a topic of informal discussion during the event, and even more so the disambiguation of concepts such as “reversibility”, “microscopic reversibility”, “detailed balance”, etc.
# 1 Markov Chains and First-Passages

Ken Sekimoto, Pedro Harunari and Charles Moslonka.<fnref n="1" /> The time evolution of numerous systems across scientific domains is well described by the mathematics of Markov chains. We will introduce the formalism, its predictions, and how to use it both analytically and numerically. The problem of first-passage times in continuous-time Markov chains (CTMCs) will be motivated and solved. Finally, we will connect the discussed points to thermodynamics and the discrete-time version.

## 1.1 Introduction

In this lecture we will introduce the main concepts and methods used in the study of Markov Chains.

In particular, the first part, which was conducted by Ken Sekimoto, will be centered around Markov chains in continuous time. Here our main goal is to set up different tools, such as the transition network (TN), master equations, and modified networks to adapt to the first-passage time problems.

In the second half-lecture (Subsections 1.7 and 1.8), conducted by Pedro Harunari, we are going to connect the different tools developed for continuous-time Markov chains to the discrete-time framework. We will first introduce those elements in the first part, and then we will present more recent developments about the first-passage time of transitions.

## 1.2 Continuous-time Markov Chains: basic notions

### 1.2.1 Notations and definitions

For basic notations, we will denote by {$a, b, c, \dots$} or {$a_1, a_2, \dots$} the discrete set of states, and $\hat{X}_t$ is the random variable representing the state $X_t$ of the system at time $t$, where $t \in [0; +\infty)$. The time evolution of $X$ is a stochastic process, and its history {$X_t, t \in \mathbb{R}^+$} is also a random variable. The sample space $\Omega$ is the set of all possible histories. We will usually denote $X_0$ by $a_0$.

We recall the Markov property for such a stochastic process:

**Definition 1** $X_t$ is a continuous-time Markovian process with respect to $t$ if the conditional probability $p(\hat{X}_{t+dt} = a | \hat{X}_{[0,t]})$ is independent of $X_s$ for all $s < t$.

The statistics of $\hat{X}_{t+dt}$ only depends on the realization of $\hat{X}_t$. It is usually said that the system forgets the past after every time step of length $dt$.

### 1.2.2 Transition rates

For two different states $a \neq b$, the conditional probability $p(\hat{X}_{t+dt} = b | \hat{X}_t = a)$ is of the order $O(dt)$ for a Markovian process. We denote by $R_{ba} > 0$ the proportionality coefficient, such that:

$$p(\hat{X}_{t+dt} = b | \hat{X}_t = a) = R_{ba}dt + O(dt^2). \qquad (1.1)$$

For several destinations {$a_1, a_2, a_3$} we have in a similar way:

$$p(\hat{X}_{t+dt} = a_i | \hat{X}_t = a) = R_{a_i a}dt + O(dt^2). \qquad (1.2)$$

A Markovian process is characterized by the set of states and the transition rates among

<footnote n="1">¹KS was the main lecturer; PH was the angel and lectured the second half; CM wrote this chapter.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 1.1: Possible transition configurations illustrating Eqs.(1.1) and (1.2)</figcaption></figure>

them. Note that, to the linear order $O(dt)$, the transition rates do not interfere with each other.

In physics, one may consider discrete problems derived from an underlying continuous process thanks to coarse-graining procedures. As an example, let us consider a random walker traveling in Europe, as in Fig.1.2. We only measure the country code with respect to time. Right after crossing a border, the walker has - for a short time period - a large probability of re-crossing the same border. Thus, if the temporal coarse-graining was not introduced, one might see multiple erratic transitions between two country codes before the walker finally moves far enough from the border.

<figure><img src="image_5.png" /><figcaption>Figure 1.2: Continuous trajectory of a random traveler in Europe, with the different countries shown in red.</figcaption></figure>

This discretization is typically non-Markovian. To obtain a Markovian trajectory in the new discrete-state continuous-time model, we need to weaken the time resolution of the trajectory, that is, introducing a time step $\Delta t$ such that faster phenomena are integrated over. More precisely, transitions $a \to b$ such that $R_{ba} \gtrsim (\Delta t)^{-1}$ should not appear in the discrete-state model. Thermodynamically, this state coarse-graining is equivalent to adding a heat bath to mask details. Descriptions with different resolutions can thus have different thermodynamics.

### 1.2.3 Transition networks

We use a network-or graph-representation for each Markov chain, in which the nodes are the states of the system, and the directed edges represent the non-zero transition rates. For the following part of the lecture, we consider ergodic transition networks i.e from any node, all the other nodes are reachable through directed edges. See Fig.1.4.
<figure><img src="image_6.png" /><figcaption>Figure 1.3: An example of ergodic transition network with four different states</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 1.4: Three examples of non-ergodic transition networks. The nodes colored in red represent potentially unreachable states.</figcaption></figure>

**Remark 1** In some other lectures of this School and Workshop, we encounter chemical reaction networks (CRN), where each node represents the state of particular constituent molecules, instead of the state of the whole system. This description will not be discussed in this chapter.

## 1.3 Simulating a trajectory: Gillespie's algorithm

Now that the basic notions of continuous-time Markov chains have been introduced, we may ask ourselves how to generate sample histories so that the statistical properties are verified. The main idea of such an algorithm is to generate a list of jumps at specific times, e.g $X_t = a \to a_i$ with $i \in \{1, \dots, n\}$.

A first but naive idea is to try a jump at every small time segment $\delta t$. This method may work but is not practical, as it is quite inefficient and could be approximate if $\delta t$ is too big.

A better-and exact - approach is the Gillespie algorithm. The idea is to generate a waiting time $\hat{T}$ between consecutive jumps. The probability of having $\hat{T} > \tau$ where $\tau > 0$ is:

$$p(\hat{T} > \tau) = \exp \left( - \sum_{i=1}^{n} R_{a_i a} \tau \right). \qquad (1.3)$$

**PROOF** The event $\hat{T} > \tau$ is equivalent to having no transitions during time intervals $(\frac{\tau}{M})k \le t < (\frac{\tau}{M})(k+1)$ for every $k \in \{0, 1, \dots, M-1\}$.

Thus, we have:

$$p(\hat{T} > \tau) = \left(1 - \sum_{i=1}^{n} R_{a_i a} \frac{\tau}{M}\right)^M \xrightarrow{M \to \infty} \exp \left( - \sum_{i=1}^{n} R_{a_i a} \tau \right). \quad \blacksquare$$
<figure><img src="image_8.png" /><figcaption>Figure 1.5: Choice protocol for the arrival state: the transition rates are “shrinked” to map the interval [0, 1].</figcaption></figure>

We then generate such a waiting time $\hat{T}$ through a uniform random variable $\hat{Y}$. We have:

$$p(\hat{T} > \tau) = p(e^{-\sum_{i=1}^{n} R_{a_i a} \hat{T}} < e^{-\sum_{i=1}^{n} R_{a_i a} \tau}) \quad (1.4)$$

$$= e^{-\sum_{i=1}^{n} R_{a_i a} \tau}. \quad (1.5)$$

Introducing $\hat{Y} := e^{-\sum_{i=1}^{n} R_{a_i a} \hat{T}}$ and $y := e^{-\sum_{i=1}^{n} R_{a_i a} \tau}$, we have:

$$p(\hat{Y} < y) = y \Rightarrow \hat{Y} \text{ is a uniform random variable on } [0, 1]. \quad (1.6)$$

After generating $\hat{Y}$ with a built-in function, we can find $\hat{T}$ such that $\hat{Y} = e^{-\sum_{i=1}^{n} R_{a_i a} \hat{T}}$.

To determine the arrival state, we notice that, in a Markovian process, the destination is determined at the last infinitesimal interval $dt$. We thus have:

$$p(\text{destination is } a_k) = \frac{R_{a_k a}}{\sum_{i=1}^{n} R_{a_i a}} \quad (1.7)$$

and the state can be decided with another uniform random variable on $[0, 1]$, see Fig.1.5.

**Remark 2** We can generalize this idea: given a 1D probability density $\rho(x)$, we can construct a random variable $\hat{X}$ that obeys $\rho(x)$. The cumulative probability up to $x \in \mathbb{R}$ is:

$$p(\hat{X} < x) = \int_{-\infty}^{x} \rho(x')dx'$$

Since this is equivalent to $p(\int_{-\infty}^{\hat{X}} \rho(x')dx' < \int_{-\infty}^{x} \rho(x')dx') = \int_{-\infty}^{x} \rho(x')dx'$, we can define $\hat{Y} := \int_{-\infty}^{\hat{X}} \rho(x')dx'$, which is a uniform random variable on $[0, 1]$, and find $\hat{X}$ by this relation.

**Exercise 1**<fnref n="2" /> Let $\rho(x, y) \in \mathbb{R}^2$ a probability density of a two-component random variable $\hat{Z} = (\hat{X}, \hat{Y})$. Describe a way of generating $\hat{Z}$ from two, mutually independent, uniform random variables on $[0, 1]$ denoted $(\hat{\xi}, \hat{\eta})$.

**Exercise 2** Consider the individual scores of a class of 134 students, with their rank sorting by increasing order (Fig.1.6). We denote by $r$ the rank and by $\psi(r) \in [0, 20]$ the score of the $r$-ranked student. We consider the score as a random variable.

<footnote n="2">²The solution to the exercises can be found at the end of this section.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 1.6: Score vs Rank of 134 students. The dotted curve is the actual data, and the smooth curve represents the continuous fitting score = $\psi(\text{rank} = r)$</figcaption></figure>

* Show that $p(\text{score} < \psi(r)) = r/134$

* Find the score probability density - over the interval [0, 20] - using $\psi(r)$ without inverting $\psi(r)$.

In Sec. 5, a python implementation of the algorithm is given for a simple chemical reaction network.

## 1.4 First-passage time problems

Consider a Markovian transition network such as the one depicted in Fig.1.3. Given an initial condition $\hat{X}_{t=0} = a_0$ and the ergodic hypothesis, the probability of the process $\hat{X}_t$ never visiting a state $a_i$ is zero. We can, therefore, define a time $\hat{T}_{FP}$ at which $\hat{X}_t$ visits $a_i$ for the first time. The random variable $\hat{T}_{FP}$ is called the first-passage time (usually abridged FPT [1]), and is a special case of stopping-time. We have in particular $p(\hat{T}_{FP} < +\infty) = 1$. Numerically, the sampling of $\hat{T}_{FP}$ can be done with a Gillespie algorithm. Moreover, we can obtain analytical results for its statistics, such as $p(\hat{T}_{FP} > \tau)$.

### 1.4.1 Master equation

Let us consider $N(\gg 1)$ copies of the transition network, starting at $\hat{X}_{t=0} = a_0$. For $t > 0$, each copy evolves independently. At a time $t$ we find $\simeq Np_t(a_i)$ copies in the state $\hat{X}_t = a_i$, with $0 \le p_t(a_i) \le 1$ and $\sum_{i=0}^n p_t(a_i) = 1$. Between $t$ and $t+dt$, the population i.e the probability of finding the system in a certain state, changes by $Np_{t+dt}(a_i) - Np_t(a_i)$. In parallel, counting all the possible transitions coming to the state $a_i$ from other states allow us to write the population influx as: $+\sum_{k(\neq i)}(R_{a_ia_k}dt)(Np_t(a_k))$. Likewise, the probability out-flux coming from the state $a_i$ to all of the other states is: $-\sum_{k(\neq i)}(R_{a_ka_i}dt)(Np_t(a_i))$. Therefore, we have the following equality between time-variation and flux:

$$Np_{t+dt}(a_i) - Np_t(a_i) = \sum_{k(\neq i)} (R_{a_ia_k}dt)(Np_t(a_k)) - \sum_{k(\neq i)} (R_{a_ka_i}dt)(Np_t(a_i)). \quad (1.8)$$

Dividing by $Ndt$ and taking the $N \to \infty$ limit, we obtain the so-called master equation:

$$\frac{dp_t(a_i)}{dt} = \sum_{k(\neq i)} R_{a_ia_k} p_t(a_k) - \sum_{k(\neq i)} R_{a_ka_i} p_t(a_i). \quad (1.9)$$

We can define the net probability flow from $a_i$ to $a_k$: $J_{a_ka_i} := -R_{a_ia_k} p_t(a_k) + R_{a_ka_i} p_t(a_i)$, also known as current, so that:

$$\frac{dp_t(a_i)}{dt} = - \sum_{k(\neq i)} J_{a_ka_i}. \quad (1.10)$$
The probability flow $J_{a_ka_i}$ can be seen as the difference of two semi-flows: the out-going flow

$$ \mathcal{J}_{a_k a_i} := R_{a_k a_i} p_t(a_i) \qquad (1.11) $$

and the in-coming flow

$$ \mathcal{J}_{a_i a_k} := R_{a_i a_k} p_t(a_k). \qquad (1.12) $$

Thus

$$ J_{a_k a_i} = \mathcal{J}_{a_k a_i} - \mathcal{J}_{a_i a_k}. \qquad (1.13) $$

The semi-flows characterize the effect of each individual possible transition, and are of particular importance when considering the transition network modifications that will be introduced in Section 1.5. For a pair of states $a_i$ and $a_k$, we say that the $a_i \leftrightarrow a_k$ transition is reciprocal if $\mathcal{J}_{a_k a_i} = \mathcal{J}_{a_i a_k}$. In this case, there is no net flow between them.

We have now switched from an individual history framework to a flow of population framework. We can now obtain a formal solution to the set of master equations. We regroup the state probabilities in a column vector: $\vec{p}_t := (p_t(a_0), \dots, p_t(a_n))^{\dagger}$. We also introduce the diagonal elements, such that: $R_{a_ia_i} := -\sum_{k(\neq i)} R_{a_k a_i}$. We can now write all of the master equations as a vector-matrix equation:

$$ \frac{d\vec{p}_t}{dt} = \mathbf{R} \vec{p}_t. \qquad (1.14) $$

The matrix $\mathbf{R}$ is called the rate matrix, and its off-diagonal elements are the transition rates: $(\mathbf{R})_{ki} = R_{a_ka_i}$. A formal solution for every $t$ is thus:

$$ \vec{p}_t = e^{\mathbf{R}t} \vec{p}_0. \qquad (1.15) $$

**Remark 3** We recall the definition of the exponential of a matrix $\mathbf{M}$:

$$ e^{\mathbf{M}} := \sum_{n=0}^{\infty} \frac{\mathbf{M}^n}{n!}. \qquad (1.16) $$

We can write the propagator–that is the path integral from an initial state to a particular later state (in this case $p(\hat{X}_t = a_k | X_0 = a_i)$–as:

$$ p(\hat{X}_t = a_k | X_0 = a_i) = (e^{\mathbf{R}t})_{a_k a_i}. \qquad (1.17) $$

### 1.4.2 First-Passage time from master equation

We can use the vector master equation to study the statistics of $\hat{T}_{FP}$, through the usage of absorbing boundary conditions (for more details, the reader may refer to Chapter XII of [2]). Figure 1.7 represents qualitatively the procedure. The master equation allows to generate individual trajectories up to a time $t$. In a way, we know the intersection (and the subsequent statistics) of the trajectories $\{a(t)\}$ with a vertical line of coordinate $t$ (Fig.1.7(a)). The first-passage time problem is, in this framework, a sort of reverse problem. We want to know the (first) intersection-time of the trajectories with a horizontal line representing a particular state ($a^*$ in Fig.1.7(b)). The main idea is to modify the transition network, by introducing particular absorbing states, meaning that all out-going transitions from them are removed. In the trajectory-space of Figure 1.7, that means that once a trajectory has reached the state of interest $a^*$, it becomes stationary (Fig.1.7(c)). We then solve the master equation (Eq.1.14) for the now-modified rate matrix. From this we can deduce the statistics of interest, such as $p(\hat{T}_{FP} < t)$ or the mean first-passage time. Below we detail this procedure applied to the network depicted in Fig.1.3.
<figure><img src="image_10.png" /><figcaption>Figure 1.7: (a) Trajectories generated by the master equation $\frac{d\vec{p}_t}{dt} = R\vec{p}_t$ representing a state variable $a$ with respect to time $t$. (b) Illustration of the first-passage time of each trajectory at the state $a^*$. (c) The same trajectories on the modified network, with an absorbing state at $a^*$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 1.8: (a): Example of a first-passage time problem at the state $a_2$ (in red) where $a_0$ (in blue) is the starting state. (b): Modified transition network, where $a_2$ is now an absorbing state (cross), with out-going transitions having been removed (light-gray lines).</figcaption></figure>
Now, we provide a step-by-step example procedure of network modification to obtain first-passage times. We consider the 4-states transition network depicted in Figs. 1.3 and 1.8(a), with a (4 × 4) rate matrix **R**. Our goal is to compute the statistics of the FPT reaching the state $a_2$ starting from $a_0$.

* We first remove the destination node (or nodes if we consider more than one state of arrival) of the FPT problem, in this case $a_2$, and replace it with an absorbing state, that is a state from which no transitions are possible. This transformation of the transition network is depicted in Fig.(1.8)(b). We denote the modified rate matrix by **R**\*, Eq.(1.19).

* We now consider the reduced state space, where all the absorbing states have been removed. In our example, the reduced state space is $(a_0, a_1, a_3)$, and the corresponding probability vector is $\vec{p}_t^* = (p_t^*(a_0), p_t^*(a_1), p_t^*(a_3))^\dagger$. The reduced master equation reads:

$$ \frac{d\vec{p}_t^*}{dt} = \mathbf{R}^* \vec{p}_t^* \qquad (1.18) $$

with

$$ \mathbf{R}^* = \begin{pmatrix} -R_{a_1 a_0} - R_{a_2 a_0} - R_{a_3 a_0} & 0 & R_{a_0 a_3} \\ R_{a_1 a_0} & -R_{a_2 a_1} & 0 \\ R_{a_3 a_0} & 0 & -R_{a_0 a_3} - R_{a_2 a_3} \end{pmatrix}. \qquad (1.19) $$

* From Eq.(1.18), we obtain the solution, given the initial condition $\vec{p}_0^* = \exp(\mathbf{R}^*t) \vec{p}_0^*$

* We can now compute the cumulative probability of the first-passage as an integral of the probability semi-flow towards the absorbing state over time:

$$ p(\hat{T}_{FP} < t) = \int_{0}^{t} [\mathcal{J}_{a_2 a_0}^{*} + \mathcal{J}_{a_2 a_1}^{*} + \mathcal{J}_{a_2 a_3}^{*}] ds \qquad (1.20) $$

$$ = \int_{0}^{t} [R_{a_2 a_0} p_{s}^{*}(a_0) + R_{a_2 a_1} p_{s}^{*}(a_1) + R_{a_2 a_3} p_{s}^{*}(a_3)] ds \qquad (1.21) $$

$$ = 1 - (p_t^*(a_0) + p_t^*(a_1) + p_t^*(a_3)) = 1 - \vec{1}^* \cdot \vec{p}_t^* \qquad (1.22) $$

with $\vec{1}^* = (1, 1, \dots, 1)$ the unit row vector over the reduced state space. While Eq.(1.23) is understood by the complementary event, $\hat{T}_{FP} \ge t$, it can also be derived from the first line using a kind of Gauss-Stokes theorem.

* The FPT probability density $\rho_{FP}(t)$ is then obtained from Eqs.(1.18) and (1.22):

$$ \rho_{FP}(t) = \frac{d}{dt} p(\hat{T}_{FP} < t) = -\vec{1}^* \cdot \mathbf{R}^* \cdot \vec{p}_t^*. \qquad (1.23) $$

From Eq.(1.23), we can compute the quantities of interest—for example, the mean first-passage time conditioned to the initial state $\vec{p}_0^* = (1, 0, 0)^T$:

$$ \begin{aligned} \langle \hat{T}_{FP} | \hat{X}_0 = a_0 \rangle &= \int_0^\infty t \rho_{FP}(t) dt \\ &= \int_0^\infty t \frac{d}{dt} p(\hat{T}_{FP} < t) dt = \int_0^\infty t \frac{d}{dt} [p(\hat{T}_{FP} < t) - 1] dt \\ &= [t(p(\hat{T}_{FP} < t) - 1)]_0^{+\infty} - \int_0^\infty (p(\hat{T}_{FP} < t) - 1) dt \\ &= 0 + \int_0^\infty \vec{1}^* \cdot \vec{p}_t^* dt \\ &= \int_0^\infty \vec{1}^* \cdot \exp(\mathbf{R}^*t) \cdot \vec{p}_0^* dt = -\vec{1}^* \cdot \mathbf{R}^{*-1} \cdot \vec{p}_0^* \end{aligned} \qquad (1.24) $$
where $\mathbf{R}^{*-1}$ denotes the inverse of the modified rate matrix, knowing that $\exp(\mathbf{R}^*t)\vec{p}_0^* \to \vec{0}^*$. It is possible to prove that $\mathbf{R}^*$ is invertible, but it is not a straightforward task; the reader may refer to Ref. [3], in particular its SM2, where it is shown that a survival matrix does not have any vanishing eigenvalue and therefore its determinant is nonzero. Note that the original matrix $\mathbf{R}$ is not invertible since 0 is an eigenvalue. Note also that $\vec{1}^* \cdot \mathbf{R}^* \neq 0$.

**Remark 4** We can derive this last result with a different approach attributed to Kramers [4]. We consider the reduced master equation, Eq.(1.18), complemented by a source term $J$ on the initial state $\vec{p}_0^*$. The equation reads:

$$ \frac{dp_t^*}{dt} = \mathbf{R}^* \vec{p}_t^* + J \vec{p}_0^*. \qquad (1.25) $$

The steady-state, denoted $\vec{p}_\infty^*$, is:

$$ \vec{p}_\infty^* = -J \mathbf{R}^{*-1} \vec{p}_0^*. \qquad (1.26) $$

By multiplying $\frac{1}{J}\vec{1}^*$ from the left,

$$ \frac{1}{J} \vec{1}^* \cdot \vec{p}_\infty^* = -\vec{1}^* \cdot \mathbf{R}^{*-1} \cdot \vec{p}_0^* = \langle \hat{T}_{FP} | \hat{X}_0 = a_0 \rangle. \qquad (1.27) $$

*Reminder : Linear algebra of master equation* We recall the spectral decomposition of the diagonalizable matrix $\mathbf{R}$: $\mathbf{R} = \mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^{-1}$ where $\mathbf{\Lambda} = \text{diag}(\lambda_1, \lambda_2, \dots, \lambda_n)$ is the diagonal matrix of eigenvalues, and $\mathbf{Q}$ is the eigenbasis representation matrix.

<figure><img src="image_12.png" /><figcaption>Figure 1.9: Columns and rows of the matrices $\mathbf{Q}$ and $\mathbf{Q}^{-1}$.</figcaption></figure>

We have $\mathbf{RQ} = \mathbf{Q\Lambda}$, and the columns of $\mathbf{Q}$ denoted by $\vec{v}_\mu$ are the right-eigenvectors of $\mathbf{R}$. Similarly, we have $\mathbf{Q}^{-1}\mathbf{R} = \mathbf{\Lambda Q}^{-1}$, so the rows of $\mathbf{Q}^{-1}$, denoted by $\vec{u}_\nu$ are the (row) left-eigenvectors of $\mathbf{R}$. We have: $\mathbf{Q}^{-1}\mathbf{Q} = I \Leftrightarrow \vec{u}_\nu \cdot \vec{v}_\mu = \delta_{\nu\mu}$ (orthonormality of the dual bases) and $\mathbf{QQ}^{-1} = \sum_\mu \vec{v}_\mu \cdot \vec{u}_\mu = I$ (completeness). The latter follows from the former. Since

$$ \begin{aligned} \exp(\mathbf{R}t) &= \exp(\mathbf{Q\Lambda Q}^{-1}t) \\ &= \sum_{n=0}^{+\infty} \frac{(\mathbf{Q\Lambda Q}^{-1})^{n}t^n}{n!} \\ &= \mathbf{Q} \left( \sum_{n=0}^{+\infty} \frac{\mathbf{\Lambda}^n t^n}{n!} \right) \mathbf{Q}^{-1} = \mathbf{Q} \exp(\mathbf{\Lambda}t) \mathbf{Q}^{-1}, \end{aligned} \qquad (1.28) $$

we have

$$ \exp(\mathbf{R}t) = \sum_{\mu} \vec{v}_{\mu} e^{\lambda_{\mu} t} \vec{u}_{\mu}. \qquad (1.29) $$

**Remark 5** Numerically the exponential $e^{\mathbf{R}^*t}$ is computed using spectral decomposition, $\mathbf{R}^* = \sum_\mu \vec{v}_\mu^* \lambda_\mu \vec{u}_\mu^*$, that is, $e^{\mathbf{R}^*t} = \sum_\mu \vec{v}_\mu^* e^{\lambda_\mu t} \vec{u}_\mu^*$. Once $e^{\mathbf{R}^*t}$ is obtained, $\vec{p}_t^*$ is given by the matrix-vector product, $\vec{p}_t^* = e^{\mathbf{R}^*t} \vec{p}_0^*$
* The matrix **R** or **R**\* can have complex eigenvalues.

* All eigenvalues of **R**\* must have strictly negative real part because $\vec{p}_t^*$ with any initial state $a_0$ should decay to the reduced zero-vector, $\vec{0}^*$, for $t \to \infty$.

* The (non-reduced) rate matrix **R** must have at least one null eigenvalue: the steady-state distribution, $\vec{\nu}_0$, satisfies $\mathbf{R}\vec{\nu}_0 = 0$. The corresponding left null eigenvector $\vec{\nu}_0$ has all components 1, a consequence of the conservation of the total probability: for all $\vec{p}_t$, we have $d/dt\vec{\nu}_0 \cdot \vec{p}_t = \vec{\nu}_0 \cdot \mathbf{R} \cdot \vec{p}_t = 0$. For further spectral properties of the rate matrix, the reader may refer to Perron-Frobenius theorems and their consequences [5].

The FPT problem shows how network modifications can be used to derive quantities of interest of the original network. We can sometimes add nodes or even replicate the original network so as to adapt to more advanced first-passage problems. Below, we study different event statistics on the same model of network, through more convoluted semi-flow re-directions.

## 1.5 Other first event problems

### 1.5.1 First and second transition time

<figure><img src="image_13.png" /><figcaption>Figure 1.10: (a): Original transition network, with the transition of interest $a_3 \to a_0$. (b): Modified TN to study the first-transition time from $a_3$ to $a_0$, with an additional absorbing state (crossed state $a_0^*$ on the right)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 1.11: Modified TN to study the second transition time from $a_3$ to $a_0$, with a replica ($a'_i$ states) and an absorbing state (crossed state $a_0^*$ on the right).</figcaption></figure>

We are now interested in the first *transition* time, which is the first time a specific transition happens in the network. In this case, we can cut the transition of inter-
est, and replace it with a transition to an absorbing state. In this case, the semi-flow,
$\mathcal{J}_{a_0 a_3}^* = R_{a_0 a_3} p_t^*(a_3')$, gives the transition time statistics.

We can extend this strategy to other transition orders, such as *second* transition times, as shown in Fig.1.10. In this case, we make use of replicas of the system without the transition of interest, linking the two states to account for the first transitions. In the example of Fig.1.10, the second transition statistics are computed with one replica, linking the states *a*₃ and the state *a*₀ of the replica, and then an absorbing state.

1.5.2 FPT with competitions

<figure><img src="image_15.png" /><figcaption>Figure 1.12: (a): Original transition network, starting from $a_0$. (b): Modified TN to study the first-passage time to $a_3$ before reaching $a_2$.</figcaption></figure>

We can also put competitions on first-passage times. For example, in Fig.1.12 we
study the first-passage time to *a*₃ before the state *a*₂ is reached. In this case we replace
*both* states by absorbing states. Then, the modified semi-flow J<sub>a₃a₀</sub>ª = R<sub>a₃a₀</sub>p<sub>t</sub><sup>*</sup>(a₀), gives
the statistics in question. We can evolve this idea to study the first-transition time for
a₃ → a₂ before the reversed transition a₂ → a₃. By placing two absorbing states at the
end of opposite transitions (Fig. 1.13), one can assess the occurrence of one or the other
(see Section 1.8).

<figure><img src="image_16.png" /><figcaption>Figure 1.13: (a) Original TN in which we study the relative first occurrence of transition between $a_2$ and $a_3$ (in light pink). (b) The modified TN where two absorbing states have been added to measure the individual semi-flows.</figcaption></figure>

1.5.3 FPT of consecutive transitions

As a last example, we can study FPT of *consecutive transitions*. For example (Fig.1.14)
we may want to study the sequence of transitions *a*₂ → *a*₃ → *a*₀. In this case, we keep
the original states *a*₃ and *a*₂, but we add a replica of the intermediate state *a*₃ and an
<figure><img src="image_17.png" /><figcaption>Figure 1.14: (a): Original transition network. (b): Modified TN to study the consecutive transitions $a_2 \rightarrow a_3 \rightarrow a_0$. In this case we compute the probability flow from $a'_3$ to $a'_0$.</figcaption></figure>

absorbing state $a'_0$. Then, the probability flow from $a'_3$ to $a'_0$ gives the FPT density of the consecutive transitions.

## 1.6 “Markovianization” of non-Markov problems

Those network modification and state replication techniques can also be used to perform exact “Markovianization” of non-Markov problems [6]. For example, we consider a 1D

<figure><img src="image_18.png" /><figcaption>Figure 1.15: Scheme of a 1D chain network on which a non-Markov jump process occurs.</figcaption></figure>

chain of states with jump probabilities depending on the actual state and the previous jump realization. We define the stochastic jump $\hat{X}_k$ that takes the value +1 when the walker goes from a state $n$ to the state $n+1$, or the value -1 if the walker goes from $n$ to $n-1$ at stage $k$. We suppose that the probability of having $(\hat{X}_k = \pm 1)$ depends on the value of the previous jump. The jump probabilities $p(\hat{X}_k = \pm 1 | X_{k-1} = \pm 1)$ are defined as follows:

<table><thead><tr><th rowspan="2">Condition</th><th colspan="2">Jump probability</th></tr><tr><th>k-th jump = +1</th><th>k-th jump = -1</th></tr></thead><tbody><tr><td>(k - 1)-th jump = +1</td><td>(1 - θ)w</td><td>θw</td></tr><tr><td>(k - 1)-th jump = -1</td><td>θw</td><td>(1 - θ)w</td></tr></tbody></table>

We can map this problem to a Markovian transition network, via a replica of the whole chain, to allow the conditions on $\hat{X}_{k+1}$ to be taken into account.

In the example depicted in Fig.1.16, the states $\{n'\}$ are only reached **from the left** to ensure that the previous transition was +1. Similarly, the states $\{n''\}$ are reached **from the right**, ensuring that the previous transition was -1. We then link the states with transition rates in accordance with the table to fully describe the Markovian network. We can then apply the previous techniques to obtain, for example, FPT statistics. Note that $n'$ and $n''$ both describe the same original state $n$, and their statistics must be added to obtain the original one.
<figure><img src="image_19.png" /><figcaption>Figure 1.16: Modified transition network derived from Fig.1.15. The non-Markovian process is rendered Markovian thanks to the double chains.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 1.17: Typical trajectories for the three different cases. Note that the smaller $\theta$ is, the more spread trajectories are. This effect is explained qualitatively by the dependency of the diffusion coefficient $D$ with $\theta$ (Eq.1.34).</figcaption></figure>

The random walk exhibits a directional persistence when the parameter $\theta$ differs from 1/2. The $0 < \theta < 1/2$ case typically represents “run and tumble” processes, whereas $1/2 < \theta < 1$ represent turn alterations (typical trajectories are depicted in Fig.1.17). Otherwise, if $\theta = 1/2$, the process is a Markovian random walk. This “fine-graining” method of Markovianization preserves the time resolution of the underlying process - as opposed to coarse-graining methods (discussed above in Sec.1.2.2).

**Remark 6** The topology of the TN in Fig.1.16 is isomorphic with that of the run-and-tumble model of an active swimmer [6]. In the latter model the nodes $n'$ and $(n - 1'')$ in Fig.1.16 is associated to the position $n$.

**Calculation of the diffusion coefficient** We take, as a time-unit, the period between consecutive tumbling. (The step displacement of the tumbling will be taken into account later as the first displacement of this period.) Since the tumbling transition rate is $\theta w$, the probability density of this period is $P(t)dt = (\theta w)e^{-\theta w t}dt$. During this period, a number $n$ of forward jumps occurs with a Poissonian distribution $P(n|t) = \frac{e^{-(1-\theta)w t}}{n!}[(1-\theta)w t]^n$. The joint probability is then,

$$P(t, n)dt = \frac{e^{-(1-\theta)w t}}{n!}[(1-\theta)w t]^n(\theta w)e^{-\theta w t}dt. \quad (1.30)$$
Note that, using the gamma function identity $\Gamma(n) = \int_0^\infty e^{-x}x^{n-1}dx = (n-1)!$ with $n \in \mathbb{N}^*$, we have:

$$ 
\begin{aligned}
P(n) &= \int_0^\infty P(t, n)dt \\
&= \int_0^\infty \frac{e^{-(1-\theta)wt}}{n!}[(1-\theta)wt]^n(\theta w)e^{-\theta wt}dt \\
&= (1-\theta)^n\theta \int_0^\infty \frac{e^{-wt}}{n!}(wt)^n wdt \\
&= (1-\theta)^n\theta \frac{1}{n!} \int_0^\infty e^x x^n dx = (1-\theta)^n\theta \quad \text{with } x = wt. \quad (1.31)
\end{aligned}
$$

The run length $\ell_n$ is given by $\ell_n = n+1$, in taking account of the initial unit step associated with the tumbling. For the $N$ pairs of periods of forward and backward runs, the total time is $T = (\sum_{k=1}^N t_k^+)+(\sum_{k=1}^N t_k^-)$ and the total displacement is $X = (\sum_{k=1}^N \ell_k^+)-(\sum_{k=1}^N \ell_k^-)$. For each period we apply $P(t_k^\pm, n_k^\pm)$, where $\ell_k^\pm = n_k^\pm + 1$. The quantity of interest is the ratio $\langle X^2 \rangle/(2\langle T \rangle)$, where $\langle \cdot \rangle$ denotes the expected value, or average. Since different $\ell_k^\pm$'s are mutually independent,

$$ \langle T \rangle = 2N \langle t \rangle \quad (1.32) $$

$$ \langle X^2 \rangle = 2N(E\langle l^2 \rangle - \langle l \rangle^2) = 2N(\langle n^2 \rangle - \langle n \rangle^2). \quad (1.33) $$

From $P(n)$ we find $\langle n^2 \rangle - \langle n \rangle^2 = \frac{(1-\theta)}{\theta^2}$, while from $P(t)$ we find $\langle t \rangle = \frac{1}{\theta\omega}$. Altogether, the diffusion constant is

$$ D = \frac{\langle X^2 \rangle}{2 \langle T \rangle} = \frac{\langle n^2 \rangle - \langle n \rangle^2}{2 \langle t \rangle} = \frac{w(1-\theta)}{2\theta}. \quad (1.34) $$

Here $w$ is the total frequency of jump, either rightward or leftward.

## 1.7 Discrete-time Markov Chains

### 1.7.1 A link between discrete and continuous Markov chains

<figure><img src="image_21.png" /><figcaption>Figure 1.18: Discretization of time on the transition network. Stationary transitions are shown in red.</figcaption></figure>

Let us consider the previously studied transition network, see Fig.1.18. The dots represent the different states of the systems and the arrows are the possible transitions. In discrete-time, the network has a similar shape, but the iteration of steps creates the possibility of sojourns when the system remains in the same state<fnref n="3" />. Those transitions

<footnote n="3">³Some transitions may happen at a finer resolution, but are not visible at the present level of description (see Section 1.2.2).</footnote>
are added to the graph, symbolized by the circular arrows. We now consider time as a discrete variable: $t = n\Delta t$ with $n \in \mathbb{N}$ and $\Delta t$ is the unitary step duration. We define the stochastic (or transfer) 4 matrix $\mathbf{P}$ by

$$ \mathbf{P} := \mathbf{I} + \Delta t \mathbf{R}, \qquad (1.35) $$

where $\mathbf{R}$ is the rate matrix for the continuous-time master equation. $\mathbf{P}$ defines how the population probability is dynamically flowing in the network. The off-diagonal elements are given by:

$$ (P)_{ab} = \Delta t(R)_{ab} = p(\hat{X}_{n+1} = a | \hat{X}_n = b) \qquad (1.36) $$

for $a$ and $b$ two states of the system. Since $(R)_{ab}$ is the probability per unit time that the system jumps from $b$ to $a$, $(P)_{ab}$ is the probability of jumping to state $a$ at time $(n+1)\Delta t$ given that state are in $b$ at time $n\Delta t$. Here $\hat{X}_n$ represents the random variable associated to the state occupied at time $n$, and $p$ is the probability measure.

The diagonal elements are:

$$ (P)_{aa} = 1 + \Delta t(R)_{aa} = 1 - \Delta t \sum_{b(\neq a)} (R)_{ba}. \qquad (1.37) $$

The diagonal element $(R)_{aa}$ is, according to the continuous-time framework, the negative sum of all transition rates from state $a$ (this results comes from the conservation of the probability norm during processes). $P_{aa}$ is thus the probability of not jumping at all:

$$ (P)_{aa} = p(\hat{X}_{n+1} = a | \hat{X}_n = a). \qquad (1.38) $$

Note that those transitions do not appear explicitly in the continuous-time framework.

The evolution equation of the probability vector $\vec{p}$ over the state space is:

$$ \vec{p}_{n+1} = \mathbf{P}\vec{p}_n \qquad (1.39) $$

which gives element-wise:

$$ p_{n+1}(b) = \sum_{a} (P)_{ba} p_n(a) \qquad (1.40) $$

$$ = \sum_{a} p(\hat{X}_{n+1} = b | \hat{X}_n = a) p_n(a) \qquad (1.41) $$

$$ = \sum_{a(\neq b)} (P)_{ba} p_n(a) + (P)_{bb} p_n(b). \qquad (1.42) $$

From Eqs. (1.36) and (1.37), we obtain:

$$ p_{n+1}(b) = \sum_{a(\neq b)} \Delta t(R)_{ba} p_n(a) + [1 + \Delta t(R)_{bb}] p_n(b) \qquad (1.43) $$

or, equivalently:

$$ \frac{p_{n+1}(b) - p_n(b)}{\Delta t} = \sum_{a(\neq b)} (R)_{ba} p_n(a) - \sum_{a(\neq b)} (R)_{ab} p_n(b). \qquad (1.44) $$

Taking the limit $\Delta t \to 0^+$ in Eq. (1.44) with $n = t$ allow us to verify the continuous-time master equation $d\vec{p}/dt = \mathbf{R}\vec{p}$. For this reason, everything said in the discrete-time

<footnote n="4">⁴The value of $\Delta t$ has to be small-enough to ensure the nonnegativity of all its entries, which is a technical detail that will not affect the connection to CTMCs.</footnote>
framework holds in the continuous-time limit, given that we take a small enough time step $\Delta t$. This equivalence is sometimes used in the other direction, switching from continuous to discrete time to prevent dealing with exponentially-distributed time intervals. This shift of description can make proofs easier with the results holding in both cases as their are a time limit away.

Now we discuss the solution of the evolution equation. In order to write a propagator for the evolution equation, we apply Eq. (1.39) n times:

$$ \begin{aligned} \vec{p}_n &= \mathbf{P}^n \vec{p}_0 \\ &= (\mathbf{I} + \Delta t \mathbf{R})^n \vec{p}_0. \end{aligned} \qquad (1.45) $$

$\mathbf{P}^n$ in the above is, therefore, the propagator. If we let $\Delta t \to 0^+$ with $n\Delta t$ fixed, we again recover the exponential propagator of the continuous-time formulation:

$$ \vec{p}_n = (\mathbf{I} + \Delta t \mathbf{R})^{t/\Delta t} \vec{p}_0 \xrightarrow[\Delta t \to 0^+]{\phantom{\Delta t \to 0^+}} e^{\mathbf{R}t} \vec{p}_0. \qquad (1.46) $$

In the following table we sum up the differences between the continuous and discrete formulations.

<table><thead><tr><th>Element</th><th>Continuous framework</th><th>Discrete framework</th></tr></thead><tbody><tr><td>Stochastic Matrix</td><td>Rate Matrix <b>R</b></td><td>Transition probability matrix <b>P</b></td></tr><tr><td>Dynamics</td><td>Master Equation <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 0.2em;">d<i>p</i><sub>t</sub></div><div style="padding: 0 0.2em;">dt</div></div> = <b>R</b><i>p</i></td><td>Evolution Equation <i>p</i><sub>n+1</sub> = <b>P</b><i>p</i><sub>n</sub></td></tr><tr><td><i>p</i>(<i>X</i><sub>t+Δt</sub> = <i>b</i>|<i>X</i><sub>t</sub> = <i>a</i>)</td><td>(<i>R</i>)<sub>ba</sub>Δ<i>t</i></td><td>(<i>P</i>)<sub>ba</sub></td></tr><tr><td>Diagonal elements</td><td>(<i>R</i>)<sub>aa</sub> = −∑<sub>b(≠a)</sub>(<i>R</i>)<sub>ba</sub> ≤ 0</td><td>0 ≤ (<i>P</i>)<sub>aa</sub> ≤ 1</td></tr><tr><td>Propagator</td><td>(<i>e</i><sup><b>R</b><i>t</i></sup>)<sub>ba</sub></td><td>(<b>P</b><sup>n</sup>)<sub>ba</sub></td></tr><tr><td>Conservation of prob.</td><td>∑<sub>b</sub>(<i>R</i>)<sub>ba</sub> = 0</td><td>∑<sub>b</sub>(<i>P</i>)<sub>ba</sub> = 1</td></tr></tbody></table>

Note that the elements of **R** have the dimension of the inverse of time whereas the elements of **P** are dimensionless. It is also worth mentioning that discrete-time processes are simpler to simulate in a computer program as we do not have to draw random time steps with a Gillespie algorithm (as described in Section 1.3). In the discrete case, $\Delta t$ is fixed and the operator only has to draw the probabilities ($P_{ba}$) from a uniform distribution between 0 and 1.

### 1.7.2 First-passage times

In a similar way as the continuous-time framework, we can compute the first-passage times by introducing absorbing states in the network - as shown in Fig.1.19. The only difference is that the stationary transition (dashed red line) cannot be removed like the other escape probabilities (Fig. 1.19b). Except for this remark, we refer to Section 1.4 for the computation of first-passage times. Note that for the problem to be properly defined, we need to assume the ergodicity of the considered network.

## 1.8 First-transition times and hidden networks

### 1.8.1 Introduction and experimental setups

In this subsection, we are no longer interested in the time when a state is reached for the first time. Instead, we want to measure the first time that a specific transition is realized in the network. This is still a very active topic of research and the reader may refer to the recent articles [7] and [8] for more details. As an illustration, consider the popular model of
<figure><img src="image_22.png" /><figcaption>Figure 1.19: In order to compute the first-passage time to the square state from the open circle state the network (a) is modified. The arrival state is replaced by an absorbing state where the escape probabilities to other states are 0 (b).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 1.20: (left) Schematic of a molecular motor walking along a microtubule, consuming ATP as an energy source and releasing ADP as a byproduct. (right) Bicyclic Kinesin model, where the states and the dashed transitions are hidden. Only the transition along the solid line is visible.</figcaption></figure>

a molecular motor moving along a microtubule and the monitoring of its activity (Fig1.20 (Left)). The only visible change in this system is the forward or backward movement of the motor. Indeed the whole metabolism behind is quite difficult experimentally to monitor precisely, and it is impossible to observe each molecule of ATP and ADP going in and out of the system. By contrast, the movement of the molecular motor itself can be measured accurately in experiments. The different states of the molecular motor are modeled by the Kinesin model. In this model, the different metabolic steps are represented by transitions in a network. Yet, the only *visible* transition is the movement of the motor, symbolized by the solid green line in Fig.1.20(right). All the other transitions and states are hidden from the observer.

Another experimental setup is a system of two quantum dots between two electrodes (Fig.1.21(a)). Each quantum dot can be occupied by one electron coming either from the electrodes or the other dot. The states of the system consist of vacant or occupied quantum dots denoted, respectively, by 0 or 1. For the two dots, the possible states are 00, 01, 10 and 11. The electrons are free to move in the system and to enter and leave through the electrodes. If we place only one detector between the two quantum dots, we do not have access to all the possible transitions (Fig.1.21(b)). In this case only the 01 ↔ 10 transition is visible.
<figure><img src="image_24.png" /><figcaption>Figure 1.21: Quantum dot experiments. Only one detector is placed (a), so that only a part of the possible transitions are observable. The equivalent network in shown on the right (b). The plain red transition is visible whereas the dotted ones are not.</figcaption></figure>

1.8.2 Inference from transition statistics

We now consider only the visible transition and its occurrences in a more formal manner. The observed trajectory (history) Γ is therefore a list of waiting times $t_i$ (with $i \in \mathbb{N}$) between occurrences of visible transition $l_i$, starting at $t = 0$. See Fig.1.22 for an illustration.

<figure><img src="image_25.png" /><figcaption>Figure 1.22: Example of a possible trajectory Γ.</figcaption></figure>

In a similar manner as in Section 1.5, we can map the question of the first occurrence of a specified transition to the question of the first occupation of a specific state, by again modifying the transition network in a particular way with absorbing states, as shown in Fig.1.23. Note that the visible transition rates are kept the same in the modified network. In this case, the first-transition time will be the first-passage time to one of the absorbing states, whatever happens in the hidden part of the system.

Considering a trajectory Γ with observed transitions {$l_i$}$_{i \in \mathbb{N}}$, we compute the first-transition time of $l_{i+1}$ given that we start (at $t = 0$) at the end of the preceding transition, $l_i$. With the modification trick in mind, we can now write the density of first-passage time according to Eq.(1.23) of Section 1.4.2:

$$ \rho_{\text{PT}}(t) = R_{l_{i+1}} p_t^*(l_{i+1} | l_i(t=0)), \qquad (1.47) $$

where $R_l$ is the transition rate of $l$. Here, $p_t^*(l_{i+1} | l_i(t=0))$ denotes the probability that - in the modified transition network - the system is in the state located at the tail of the transition $l_{i+1}$ at time $t$, given that the system is in the state at the tip (head) of transition $l_i$ at time $t = 0$. The system can, therefore, perform any succession of hidden transition during 0 and $t$ with no visible transitions occurring, because the absorbing states prevent those events. The time $t$ is hence called the “intertransition time”.

From the statistics of the intertransition time (which are experimentally measurable quantities), one can probe characteristics of the whole system, such as the amount of
<figure><img src="image_26.png" /><figcaption>Figure 1.23: Example of a transition network with hidden states (dotted filled circles) and transitions (dotted curves), and one visible transition (solid curves). The network can be modified to switch from a first-transition time problem to a first-passage time problem, thanks to the addition of absorbing states. Note that, in this case, these absorbing states cannot be physically-meaningful initial positions, as they are only computational artifacts.</figcaption></figure>

disorder of the molecular system, bounds for the entropy production rate and the thermo-dynamic efficiency, insights of the hidden network topology and possibly more quantities that are yet to be found [7, 8].
## Solution to the exercises

**Exercise 1** Let us first take care of the $x$ component. Similar to the 1D case, we define $\hat{X}$ such that:

$$ \hat{\xi} = \int_{-\infty}^{\hat{X}} \rho_x(x) dx \qquad (1.48) $$

where $\rho_x(x)$ is the marginal $\rho_x(x) = \int_{\mathbb{R}} \rho(x, y) dy$. In practice, one needs the inverse function of the cumulative probability distribution which can be computed analytically or numerically.

For a given $\hat{X}$, we define $\hat{Y}$ such that:

$$ \hat{\eta} = \int_{-\infty}^{\hat{Y}} \rho(y|\hat{X}) dy \qquad (1.49) $$

where $\rho(y|x)$ denotes the conditional probability density: $\rho(y|x) = \rho(x, y)/\rho_x(x)$. Then, for given $(x, y)$, the probability attributed to the surface element $(dx \wedge dy)$ by the random variables $\xi$ and $\eta$ is:

$$ d\xi \wedge d\eta = \rho_x(x)dx \wedge \left( \rho(y|x)dy + \left[ \int_{-\infty}^{y} \frac{\partial}{\partial x} \rho(y'|x)dy \right] dx \right) \qquad (1.50) $$

$$ = \rho_x(x)\rho(y|x)(dx \wedge dy) + 0 = \rho(x, y)(dx \wedge dy) \qquad (1.51) $$

because $(dx \wedge dx) = 0$. Therefore, $\hat{X}$ and $\hat{Y}$ corresponds to the density $\rho(x, y)$.

## Exercise 2

* Question 1: This is the definition of a ranked series. Since the interval of ranks is [0, 134], the fraction $r/134$ is the number of scores below the score of the $r$-th student, which is $\psi(r)$. Hence $r/134 = p(\hat{S} < \psi(r))$, where $\hat{S}$ is the random variable corresponding to the score of a randomly chosen student.

* Question 2: By taking the derivative of the relation, $p(\text{score} < \psi(r)) = r/134$, with respect to $r$, we have $\rho(\psi(r))\psi'(r) = (134)^{-1}$. Then the parametric presentation, $(134\psi'(r))^{-1}$ vs $\psi(r)$ gives the relation $\rho(s)$ vs $s$.
## 2 Network Thermodynamics

Sara Dal Cengio, Nahuel Freitas and Paul Raux.<fnref n="5" /> We introduce the basic concepts and computational tools of network theory. We discuss their use in the thermodynamics of Markov jump processes, including the Markov chain tree theorem and the Schnakenberg decomposition of forces and fluxes. We draw connections with the deterministic description of electrical circuits and chemical reaction networks, highlighting similarities and differences.

### 2.1 Introduction

This chapter is dedicated to motivate the use of graph theory in the field of nonequilibrium thermodynamics. Graph theory found its first application in electricity. Therein, Kirchhoff gave a major contribution to the field by stating the current and the voltage laws, opening up to the understanding of electrical circuits as nonequilibrium systems involving one (type of) potential and current. Since then, graph theory has expanded its scope. Generalised Kirchhoff Current Law and Kirchhoff Voltage Law have been developed to accommodate systems with multiple currents and potentials, such as chemical networks. A vast recent literature on this subject can be found [9–14], but the review by Schnakenberg in 1976 [15] is considered as seminal in the field. In this lecture, we focus on two representatives of this literature: the Markov chain tree theorem and the decomposition of observables on a graph. These notes are organised as follows: in the first subsection, we introduce the elements of graph theory needed to describe a Markov chain on a graph; in the second subsection we recall some elements of linear algebra and use them to prove the Markov chain tree theorem; the last subsection is dedicated to the decomposition of physical observables on the graph.

### 2.2 Elements of graph theory

In this first subsection, we introduce the few basic elements of graph theory that we need in order to illustrate the relevance of graph theory in nonequilibrium thermodynamics. Across this subsection, we will illustrate the definitions on the four vertex graph already used as an example in [15] (See Figure 2.1).

#### 2.2.1 Graph, incidence matrix and spanning tree

Let $G(N, E)$ be a connected graph with $N$ vertices and $E$ edges. We denote respectively $V$ and $E$, the vertex space and the edge space of $G$, so that $\dim(V) = N$ and $\dim(E) = E$. We assign to each edge $e$ an arbitrary orientation by defining a source vertex $s(e)$ and a target vertex $t(e)$, so that an edge can be more precisely denoted as $e \equiv (s(e) \to t(e))$. We emphasize that this orientation is a convention used to describe the network and does not convey any physical meaning concerning the directionality of transitions. Indeed, in the physical cases, transitions along an edge $e$ can be performed in both directions according to micro-reversibility. We denote $-e$ the reverse direction: $-e \equiv (t(e) \to s(e))$.

**Incidence matrix** The incidence matrix $D$, characterizing $G$, is of size $N \times E$, and is defined as follows:

$$D_{i,e} = \delta_{i,t(e)} - \delta_{i,s(e)}, \qquad (2.1)$$

<footnote n="5">⁵SDC was the lecturer; NF was the angel; PR wrote this chapter.</footnote>
where $\delta_{i,j}$ denotes the Kronecker delta, and the indices $i$ and $e$ span respectively the set of vertices and edges of $\mathcal{G}$. Eq. (2.1) encodes the fact that each edge has exactly one source and one target. As a consequence, the row vector $\ell^\top = (1, \dots, 1)$ is a left null vector of $D$, so that:

$$ \sum_i \ell_i D_{ie} = \sum_i D_{ie} = 0, \forall e. \qquad (2.2) $$

It follows that the lines of $D$ are not linearly independent and the rank of the incidence matrix, rank$D$, is strictly less than $N$. Specifically, since graph $\mathcal{G}$ is connected,

$$ \mathrm{rank} D = N - 1 \qquad (2.3) $$

**Proof:**

* First, the rank nullity theorem of linear algebra relates rank$D$ with the dimension of the kernel, $\dim(\ker D)$, so that:

$$ \mathrm{rank} D = \mathrm{rank} D^\top = N - \dim(\ker D^\top) \qquad (2.4) $$

where $D^\top$ is the transpose of matrix $D$. Thus, to prove that $\mathrm{rank} D = N - 1$ is fully equivalent to show that $\dim(\ker D^\top) = 1$.

* We have seen already that the vector $\ell^\top$ belongs to $\ker D^\top$. Ad absurdum, let assume that there exists another vector, independent from $\ell$, which belongs to $\ker D^\top$. If such a vector exists, one could build, by linear combination with $\ell$, a non-zero vector $\ell' \in \ker D^\top$ containing a 0 component $\ell'_i = 0$; this is impossible since, using Eq. (2.2) by recursion along the connected graph $\mathcal{G}$, starting from vertex $i$, we find $\ell'_j = 0 \forall j$. This is against the original assumption $\ell' \neq 0$; thus the kernel of $D$ is uniquely spanned by $\ell$ and $\mathrm{rank} D = N - 1$.

<figure><img src="image_27.png" /><figcaption>Figure 2.1: This graph has four vertices and five edges. The vertex space is $\mathcal{V} = (1, 2, 3, 4)$. The edge space is $\mathcal{E} = (e_1, e_2, e_3, e_4, e_5)$. Consider for example edge $e_1$, then the arbitrary orientation chosen here gives: $s(e_1) = 2$ and $t(e_1) = 1$ so that $e_1$ is the edge connecting $2 \to 1$.</figcaption></figure>

**Example:** Figure 2.1 shows an oriented graph with four vertices and five edges. From now on, we will denote it $\mathcal{G}_{\mathrm{ex}}$. The vertex space and the edge space are respectively given by $\mathcal{V} = (1, 2, 3, 4)$ and $\mathcal{E} = (e_1, e_2, e_3, e_4, e_5)$. With this ordering of the vertices and the
edges, the incidence matrix of $\mathcal{G}_{ex}$ reads:

$$D = \begin{pmatrix} 1 & 0 & 0 & 0 & 1 \\ -1 & 1 & 1 & 0 & 0 \\ 0 & -1 & 0 & 1 & -1 \\ 0 & 0 & -1 & -1 & 0 \end{pmatrix} \qquad (2.5)$$

The first column in $D$ corresponds to $e_1$. It has only two non-zero entries:

* a 1 for the target vertex of $e_1$: $t(e_1) = 1$

* a $-1$ for the source vertex of $e_1$: $s(e_1) = 2$

**Spanning tree** Let define a spanning tree $\mathcal{T}_G$ as a sub-graph of $\mathcal{G}$ such that:

* $\mathcal{T}_G$ contains all the vertices of $\mathcal{G}$.

* All the edges in $\mathcal{T}_G$ are edges of $\mathcal{G}$.

* $\mathcal{T}_G$ is connected, i.e. there exists a succession of edges from any vertex to any other.

* $\mathcal{T}_G$ contains no loop, i.e. no cyclic sequences of edges.

Several remarks follow:

* $\mathcal{T}_G$ is not unique.

* $\mathcal{T}_G$ contains exactly $N-1 = \text{rank}D$ edges, and its corresponding incidence matrix is full column rank. This is due to the absence of loops.

* The definition of spanning tree does not specify its edge orientation. If the edges of a spanning tree are oriented all towards vertex $i$, the spanning tree is said to be rooted in $i$. We denote the set of $i$-rooted spanning trees as $\mathcal{T}^i$. In the following, we denote with $\mathcal{T}$ the set of un-directed spanning tree, i.e. spanning trees with no edge orientation (see Figure 2.2).

**Example:** The 8 un-directed spanning trees of $\mathcal{G}_{ex}$ are given in Figure 2.2. The set of 1-rooted spanning trees is given in Figure 2.3. Consider in particular the spanning tree circled in grey in the upper left corner of Figure 2.3, that we denote $\mathcal{T}_{\mathcal{G}_{ex}}^1$. $\mathcal{T}_{\mathcal{G}_{ex}}^1$ is one element of the set $\mathcal{T}^1$ and it consists of the three first edges of $\mathcal{G}$: $(e_1, e_2, e_3)$. Accordingly, the incidence matrix of this spanning tree, denoted $D$, corresponds to the three first columns of $D$ in Eq. (2.5) so that:

$$D = \begin{pmatrix} 1 & 0 & 0 \\ -1 & 1 & 1 \\ 0 & -1 & 0 \\ 0 & 0 & -1 \end{pmatrix} \qquad (2.6)$$

### 2.2.2 Cycles and cocycles

From a given choice of spanning tree $\mathcal{T}_G$, two subsets of $\mathcal{E}$ can be identified:

* The set of chords is the set of edges that do not belong to $\mathcal{T}_G$. Since removing the chords is equivalent to opening all the loops of $\mathcal{G}$, the cardinality of chords is equal to the number of loops of $\mathcal{G}$, i.e. $c = E - N + 1$. From the viewpoint of the incidence matrix, the chords correspond to $c$ dependent columns of $D$. We denote by $\alpha$ the index spanning the set of chords. In the example in Figure 2.3, the set of chords whose removal generates $\mathcal{T}_{\mathcal{G}_{ex}}^1$ is $(e_4, e_5)$.
<figure><img src="image_28.png" /><figcaption>Figure 2.2: $\mathcal{T}$: Set of un-directed spanning trees of $\mathcal{G}_{\text{ex}}$.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Figure 2.3: $\mathcal{T}^1$: Set of 1-rooted spanning trees of $\mathcal{G}_{\text{ex}}$.</figcaption></figure>

* The set of cochords is the set of edges that belong to $\mathcal{T}_G$. From the viewpoint of the incidence matrix, the cochords correspond to linearly independent columns of $\mathcal{D}$. This encodes the absence of loops in $\mathcal{T}_G$. We denote by $\gamma$ the index that spans the set of cochords. In the example of Figure. 2.3, the set of cochords of $\mathcal{T}_{\mathcal{G}_{\text{ex}}}^1$ is $(e_1, e_2, e_3)$.

In what follows, we always order the edges of $\mathcal{G}$ such that $1 \le \gamma \le N-1$ and $N \le \alpha \le E$.

**Cycle** A cycle is the loop obtained by placing back a chord in $\mathcal{T}_G$ and removing the cochords that do not belong to the loop just created. By convention, we assign to the cycle the orientation of its generating chord. Then, formally we define cycle vectors $c^\alpha \in \mathcal{E}$
<figure><img src="image_30.png" /><figcaption>Figure 2.4: Cycles of $\mathcal{G}_{ex}$ built out of $\mathcal{T}_{\mathcal{G}_{ex}}^1$.</figcaption></figure>

as vector with components:

$$ c_e^\alpha = \begin{cases} 1 & \text{if } e \text{ belongs to the cycle generated by the chord } \alpha \\ -1 & \text{if } -e \text{ belongs to the cycle generated by the chord } \alpha \\ 0 & \text{otherwise} \end{cases} \qquad (2.7) $$

Given the one-to-one correspondence between chords and cycles, we index them with the same index $\alpha$.

**Example:** Figure 2.4 shows how to obtain the cycles of $\mathcal{G}_{ex}$ out of $\mathcal{T}_{\mathcal{G}_{ex}}^1$. For instance, placing back the chord $e_4$, we obtain the cycle $c^{\alpha=4}$. The cochord $e_1$ (dashed line in Figure 2.4) does not belong to $c^{\alpha=4}$, it is thus removed. The orientation of $c^{\alpha=4}$ is given by the orientation of $e_4$ in $\mathcal{G}_{ex}$. The edges that belong to $c^{\alpha=4}$ are $(e_2, -e_3, e_4)$. Thus the cycle vector of $c^{\alpha=4}$ reads:

$$ c^{\alpha=4} = \begin{pmatrix} 0 \\ 1 \\ -1 \\ 1 \\ 0 \end{pmatrix} \qquad (2.8) $$

**Cocycle** Removing a cochord from $\mathcal{T}_G$ creates a cut, i.e. it separates $\mathcal{T}_G$ into two disconnected components called islands. We name respectively source (target) island the island that contains the source (target) vertex of the generating cochord. Then, a cocycle is defined as the set of edges in $\mathcal{G}$ which reconnects the two islands, from source to target. Accordingly the cocycle vectors $c^\gamma \in \mathcal{E}$ are formally defined as vector with components:

$$ c_e^\gamma = \begin{cases} 1 & \text{if } e \text{ belongs to the cocycle generated by the cochord } \gamma \\ -1 & \text{if } -e \text{ belongs to the cocycle generated by the cochord } \gamma \\ 0 & \text{otherwise} \end{cases} \qquad (2.9) $$
<figure><img src="image_31.png" /><figcaption>Figure 2.5: Blue edges: cocycle built out of the removal of one cochord from $\mathcal{T}_{G_{ex}}^1$.</figcaption></figure>

Given the one-to-one correspondence between cochords and cocycles, we index them with the same index $\gamma$.

**Example:** Figure 2.5 shows how to obtain a cocycle of $\mathcal{G}_{ex}$. We start by removing the cochord $e_2$ from $\mathcal{T}_{\mathcal{G}_{ex}}^1$. This creates a cut in the graph, with the two vertices of $e_2$, here $s(e_2) = 3$ and $t(e_2) = 2$, belonging to two disconnected components. Accordingly, we name source island the disconnected component containing vertex 3 and target island the component containing vertex 2. Then, the cocycle generated by $e_2$, denoted as $c^{\gamma=2}$, is the set of edges that reconnects $\mathcal{G}_{ex}$ from source to target. In this case, these edges are $(e_2, -e_4, e_5)$ and the cocycle vector of $c^{\gamma=2}$ reads:

$$c^{\gamma=2} = \begin{pmatrix} 0 \\ 1 \\ 0 \\ -1 \\ 1 \end{pmatrix} \qquad (2.10)$$

Several remarks go as follows:

* $\forall \alpha, c^{\alpha} \in \ker D$. This means that for all $\alpha$, the following property holds:

$$\sum_{e} D_{ie} c_{e}^{\alpha} = 0 \quad \forall i \qquad (2.11)$$

which encodes the fact that, along a cycle, every vertex has exactly one incoming and one outgoing edge.

* All cocycle vectors are linearly independent since, by construction, every cochord belongs to exactly one cocycle. Hence:

$$c^{\gamma} \cdot e^{\gamma'} = \delta_{\gamma, \gamma'} \quad \forall 1 \le \gamma, \gamma' \le N-1 \qquad (2.12)$$

where $e^{\gamma}$ is the canonical vector of $\mathcal{E}$ associated with the cochord $\gamma$ so that $e_e^{\gamma} = \sum_e \delta_{\gamma, e}$.
* All cycle vectors are linearly independent since, by construction, every chord belongs to exactly one cycle. Reformulating this in mathematical terms, we have:

$$ \mathbf{c}^{\alpha} \cdot \mathbf{e}^{\alpha'} = \delta_{\alpha, \alpha'} \quad \forall N \le \alpha, \alpha' \le E \qquad (2.13) $$

where $\mathbf{e}^{\alpha}$ is the canonical vector of $\mathcal{E}$ associated with the chord $\alpha$ so that $e_e^{\alpha} = \sum_e \delta_{\alpha, e}$.

* Cycles and cocycles span orthogonal sets, namely:

$$ \mathbf{c}^{\alpha} \cdot \mathbf{c}^{\gamma} = 0 \quad \forall \alpha, \gamma \qquad (2.14) $$

This follows from the definition of cycles and cocycles which ensures that $(\pm e_{\gamma}, e_{\alpha}) \in \mathbf{c}^{\alpha}$ iff $(e_{\gamma}, \mp e_{\alpha}) \in \mathbf{c}^{\gamma}$.

It follows that cycles and cocycles $(\mathbf{c}^{\gamma}, \mathbf{c}^{\alpha})$ form a basis for $\mathcal{E}$. Specifically cycles span $\ker D$ and cocycles span its orthogonal complement space, the coimage of $D$, $\mathrm{im} D^{\top}$. We hereafter summarize the major algebraic properties encountered so-far:

$$ \ker D \perp \mathrm{im} D^{\top} \qquad (2.15) $$

$$ \ker D = \mathrm{span} \mathbf{c}^{\alpha} \qquad (2.16) $$

$$ \mathrm{im} D^{\top} = \mathrm{span} \mathbf{c}^{\gamma} \qquad (2.17) $$

## 2.3 Network theory applied to CTMC: the Markov chain tree theorem

In this subsection, we are interested in the case in which the graph $\mathcal{G}$ describes a CTMC. Within this framework, it is often of interest to look for the steady-state distribution of the CTMC as it gives a first insight into the dynamics of the system. The Markov chain tree theorem [16] gives a nice interpretation of the stationary state of a Markov chain in terms of algebraic graph theory. The outline of this subsection is the following: first we recall properties of the Markovian dynamics. Then, we recollect the necessary linear algebra basics needed to prove the theorem. Lastly, we provide a proof.

### 2.3.1 Markovian dynamics in a nutshell

For a more complete introduction to Markovian dynamics and its application to stochastic thermodynamics, the reader can refer to [17].

A Markovian dynamics on a connected graph $\mathcal{G}$ satisfies the master equation:

$$ \begin{aligned} \partial_t p_i(t) &= \sum_j [r(i|j)p_j(t) - r(j|i)p_i(t)] \\ &= \sum_e D_{ie} p_e(t) \end{aligned} \qquad (2.18) $$

where $D$ is the incidence matrix of $\mathcal{G}$, $r(i|j)$ is the transition probability rate to jump from $j$ to $i$, and we have introduced the net probability flux along the oriented edge $e$ as:

$$ j_e \equiv r(e)p_{s(e)} - r(-e)p_{t(e)} \qquad (2.19) $$

where $r(e) \equiv r(t(e)|s(e))$ and $r(-e) \equiv r(s(e)|t(e))$. The second equality in Eq. (2.18) shows that, apart from a minus sign, the master equation can be interpreted as a continuity equation for the probabilities $p_i$'s, with $D$ a discrete divergence operator. The master equation can also be written as follows:

$$ \partial_t p_i(t) = \sum_j R_{ij} p_j(t) \qquad (2.20) $$
where we have introduced the $N \times N$ rate matrix defined as:

$$ \mathbf{R}_{ij} = r(i|j) - \left( \sum_k r(k|i) \right) \delta_{ij} \qquad (2.21) $$

Then, probability conservation, namely $\sum_i p_i(t) = 1 \forall t$, is ensured by the fact that the elements of the columns of $\mathbf{R}$ sum to 0:

$$ \sum_i \mathbf{R}_{ij} = 0 \forall j \qquad (2.22) $$

Eq. (2.22) is a null linear combination of the lines of $\mathbf{R}$. It implies that the row vector $\mathbf{\ell}^\top$ [that we already met] is a left null eigenvector of $\mathbf{R}$, and that $\mathbf{R}$ has determinant 0. Since the graph is connected, Perron-Frobenius theorem ensures the uniqueness of the left null eigenvector, as well as the existence of a unique stationary state [18], since $\dim(\ker \mathbf{R}) = \dim(\ker \mathbf{R}^\top) = 1$. Accordingly, this stationary solution $\mathbf{p}_\infty$ is obtained as a right null eigenvector of $\mathbf{R}$:

$$ \sum_j \mathbf{R}_{ij} p_j(\infty) = 0 \qquad (2.23) $$

Then one way to obtain the stationary probability distribution of our CTMC is to directly compute the right null eigenvector of $\mathbf{R}$. However, we will see that a graph theoretical expression of $\mathbf{p}(\infty)$ can also be found.

### 2.3.2 Markov chain tree theorem

As introduced in 2.2, let $\mathcal{T}$ be the set of un-rooted spanning tree of $\mathcal{G}$ and $\mathcal{T}^i$ be the corresponding set of $i$-rooted spanning trees. Then, we define the weight vector $\mathbf{w} \in \mathcal{V}$ with components:

$$ w_i = \sum_{\mathcal{T}} w(\mathcal{T}^i) = \sum_{\mathcal{T}} \prod_{e \in \mathcal{T}^i} r(e) \quad \forall i \qquad (2.24) $$

given by the sum over all spanning trees of the product of their transition rates with the edges oriented toward $i$. The Markov chain tree theorem states that the stationary distribution is given by the following formula:

$$ p_i(\infty) = \frac{w_i}{\sum_j w_j} \qquad (2.25) $$

**Example:** Let us illustrate the theorem by computing the stationary probability of a CTMC on a graph $\mathcal{G}$ with vertex space $\mathcal{V} = (1, 2, 3)$ and edge space $\mathcal{E} = (e_1, e_2, e_3)$ (See Figure 2.6). According to the Markov tree chain theorem, to compute the stationary probability to be in state $i$, we need to compute the weight $w_i$. To do so, we need to find all the spanning trees of $\mathcal{G}$ rooted in state $i$. Figure 2.6 shows for instance the three spanning trees rooted in state 2. Then the weight $w_2$ is obtained by summing, among the possible rooted spanning trees, the products of the transition rates. Here, for state 2, we thus have:

$$ w_2 = r(1|3)r(2|1) + r(2|1)r(2|3) + r(3|1)r(2|3) \qquad (2.26) $$

The process can be repeated for states 1 and 3 to obtain $w_1, w_3$. By normalizing Eq. (2.26), we then obtain the stationary probability to be in state 2:

$$ p_2(\infty) = \frac{r(1|3)r(2|1) + r(2|1)r(2|3) + r(3|1)r(2|1)}{Z} \qquad (2.27) $$

where $Z = \sum_i w_i$.
<figure><img src="image_32.png" /><figcaption>Figure 2.6: a) Oriented graph G with 3 vertices. b) Set of 2-rooted spanning trees of G.</figcaption></figure>

### 2.3.3 Proof of the theorem

There exists no trivial proof of the above theorem. In this subsection, we propose a compact proof for CTMC following the spirit of [19].

**Linear algebra** We start this proof subsection with a short reminder of linear algebra. Given a square matrix $A$, the cofactor matrix of $A$ denoted $C_A$ is defined as follows:

$$ (C_A)_{ij} = (-1)^{i+j} (M_A)_{ij} \qquad (2.28) $$

where $(M_A)_{ij}$ is the determinant of the matrix obtained by removing line $i$ and column $j$ from matrix $A$, also called minor of $A$. It can thus be written as $(M_A)_{ij} = \det A_{(i,j)}$. The adjugate matrix is defined as the transpose of $C_A$:

$$ \mathrm{adj} A = (C_A)^{\top} \qquad (2.29) $$

We will use the following property of the adjugate matrix:

$$ A \mathrm{adj} A = \mathrm{adj} A A = \det A \mathbf{1} \qquad (2.30) $$

where $\mathbf{1}$ is the identity matrix. Finally, we introduce the Cauchy-Binet formula that will be used later on. Let $A$ and $B$ be two conformable matrices of respective size $m \times n$ and $n \times m$. The Cauchy-Binet formula reads:

$$ \det AB = \sum_{S_m^n} \det A^{(S_m^n)} \det B^{\top(S_m^n)} \qquad (2.31) $$

where $S_m^n$ denotes the sub-ensembles of $m$ elements out of $(1, ..., n)$. The cardinality of such sub-ensembles is the binomial coefficient: $\binom{n}{m}$. Then, $A^{(S_m^n)}$ corresponds to the $m \times m$ matrix whose columns are the columns of $A$ indexed by $S_m^n$, and the same holds for $B^\top$.

**Usefulness of the adjugate matrix:** As a first lemma, we show that the diagonal of the adjugate matrix of $R$ is related to $p(\infty)$. Since $\det R = 0$, Eq. (2.30) applied to the rate matrix gives:

$$ \mathrm{R adj} R = \mathbb{0}_N \qquad (2.32) $$
where we have denoted $\mathbb{0}_N$ the null matrix of size $N \times N$. Therefore, every column of the adjugate matrix is a right null eigenvector of $R$. Given that $\dim(\ker R) = 1$, it means that all columns of $\text{adj}R$ are proportional to each other. Moreover, using the second equality of Eq. (2.30), we have:

$$R^{\top} (\text{adj}R)^{\top} = \mathbb{0}_N \qquad (2.33)$$

meaning that the column vectors of $(\text{adj}R)^{\top}$ are in $\ker R^{\top}$. This in turn implies that every vector of $\ker R^{\top}$ is proportional to $\ell$. Therefore, denoting $v^i$ the $i$-th column vector of $(\text{adj}R)^{\top}$, we have $\forall i \ v^i = \alpha_i \ell$, with $\alpha_i \in \mathbb{R}$. It follows that every columns of $\text{adj}R$ is equal to $(\alpha_1, ..., \alpha_N)^T$. In particular, they are all equal to the vector given by the diagonal entries of $\text{adj}R$. Since this vector lives in $\ker R$, we obtain from it the stationary probability by imposing normalization:

$$p_i(\infty) = \frac{(\text{adj}R)_{ii}}{\sum_i (\text{adj}R)_{ii}} \qquad (2.34)$$

We have shown how the adjugate matrix relates to the stationary distribution of our CTMC. We now turn to the task of drawing a connection to graph theory.

**Connection with graph theory:** We have shown in the last paragraph that the columns of $\text{adj}R$ are all equal. We may thus write:

$$(\text{adj}R)_{ij} = (\text{adj}R)_{ii} = (M_R)_{ii} \qquad (2.35)$$

where the second equality originates from the definition of the adjugate matrix Eq. (2.29). The connection to graph theory goes as follows. First, we have seen in subsection 2.3.1 that the dynamics of the CTMC can be described either in the space of probabilities $p_i$ with the rate matrix $R$ or in the space of probability fluxes $j_e$ with the matrix $D$. We now show how to relate the two above matrices. First we introduce the matrix $\mathcal{R}$ defined as:

$$\mathcal{R}_{i,e} = r(-e)\delta_{i,t(e)} - r(e)\delta_{i,s(e)} \qquad (2.36)$$

Then, we realize that $R$ can be written in terms of $\mathcal{R}$ and $D$ as follows:

$$R = D R^{\top} \qquad (2.37)$$

Eq. (2.37) makes explicit the relation between the rate matrix of the CTMC and the incidence matrix $D$ containing information on the topology of the underlying network. Taking back equation (2.35), we have:

$$\begin{aligned} (\text{adj}R)_{ii} &= (\text{adj}R)_{ij} = (M_R)_{ii} \\ &= \det D_{\backslash(i, \cdot)} R_{\backslash(\cdot, i)}^{\top} \\ &= \sum_{S_{N-1}^{E}} \det D_{\backslash(i, \cdot)}^{(S_{N-1}^{E})} \det R_{\backslash(\cdot, i)}^{(S_{N-1}^{E})} \\ &= \sum_{\mathcal{T}} \det D_{\backslash(i, \cdot)}^{(\mathcal{T})} \det R_{\backslash(\cdot, i)}^{(\mathcal{T})} \end{aligned} \qquad (2.38)$$

In the second equality we have replaced $R$ by its definition in terms of $D$ and $\mathcal{R}$. To compute $(M_R)_{ii}$, we have to remove line $i$ and column $i$ of $R$, which is equivalent to removing line $i$ of $D$ and column $i$ of $R^{\top}$, indicated respectively as $D_{\backslash(i, \cdot)}$ and $R_{\backslash(\cdot, i)}^{\top}$. In the third line, we have used the Cauchy-Binet formula introduced earlier. Here the sub-ensembles are the choices of $N-1$ columns of $D_{\backslash(i, \cdot)}$ and $R_{\backslash(\cdot, i)}^{\top}$. One then has to
compute the determinant of the resulting $(N-1) \times (N-1)$ square matrices. Crucially, the only choices that contribute to the determinant are the ones in which the chosen $N-1$ elements are linearly independent. In graph theory, the set of choices of independent edges is the set of spanning trees. Therefore, the sum on $S_{N-1}^E$ can be changed on a sum over $\mathcal{T}$. For all $\mathcal{T}$, matrices $D_{(i,..)}^{(\mathcal{T})}$ and $\mathcal{R}_{(.,i)}^{(\mathcal{T})}$ are full-rank and can be made upper triangular using elementary row and column operations. Notably, given the common structure of $D$ and $\mathcal{R}$ (see Eq. (2.1) and Eq. (2.36)), the required elementary operations are the same for the two matrices, leaving the sign of the product of determinants unchanged. Once the two matrices are reduced in upper triangular form, their determinant is given by the product of the diagonal entries. By an adequate labeling of the lines and columns, $D_{(i,..)}^{(\mathcal{T})}$ and $\mathcal{R}_{(.,i)}^{(\mathcal{T})}$ can be made upper triangular with respectively $-1$ and $-r(e)$ diagonal entries for every edge directed inward vertex $i$, and $+1$ and $+r(-e)$ diagonal entries for every edge directed outward vertex $i$ (see Appendix A in [20] for details). It follows that the resulting product of rates corresponds to a sum over spanning trees rooted in $i$ and we may write:

$$ (\mathrm{adj} \mathcal{R})_{ii} = \sum_{\mathcal{T}} \prod_{e \in \mathcal{T}^i} r(e) = w_i \qquad (2.39) $$

Finally, ensuring the normalization of the probabilities $p_i$'s completes the proof.

## 2.4 Network theory applied to observable decomposition

In this subsection, we apply the notions of cycles and cocycles introduced in subsubsection 2.2.2 to the decomposition of physical observables on $\mathcal{G}$. Thermodynamics imposes strong conditions on the decomposition. We show in particular that we are able to locate thermodynamic processes on subsets of $\mathcal{G}$.

### 2.4.1 Observables decomposition

**Fluxes** We first recall that $(c^\gamma, c^\alpha)$, as well as $(e^\gamma, c^\alpha)$ and $(c^\gamma, e^\alpha)$ form bases of the edge space $\mathcal{E}$, due to the orthogonality conditions in Eqs. (2.13-2.12). Therefore, since the vector of probability fluxes $j(t)$, whose components are given by Eq. (2.19), belongs to $\mathcal{E}$, we may decompose it as follows:

$$ j(t) = \sum_{\gamma=1}^{N-1} \mathcal{J}_{\gamma}(t) e^{\gamma} + \sum_{\alpha=N}^{E} \mathcal{J}_{\alpha}(t) c^{\alpha} \qquad (2.40) $$

where coefficients $\mathcal{J}_\gamma(t) \equiv j(t) \cdot c^\gamma$ and $\mathcal{J}_\alpha^c(t) \equiv j(t) \cdot e^\alpha$, and we used the orthogonality properties Eqs. (2.13-2.12). Therefore, coefficients $\mathcal{J}_\gamma(t)$ for $1 \le \gamma \le N-1$ in Eq. (2.40) correspond to the $N-1$ probability fluxes flowing across the cocycles. More can be said about these coefficients. Indeed, at steady-state, the master equation (2.18) reduces to:

$$ \mathcal{D}j(\infty) = 0 \qquad (2.41) $$

Thus, at steady state $j(\infty)$ lives in the kernel of $\mathcal{D}$. Since cycles form a basis of $\ker \mathcal{D}$, the only term that survives at steady state is the cycle term of the decomposition, so that:

$$ j(\infty) = \sum_{\alpha=N}^{E} \mathcal{J}_{\alpha}(\infty) c^{\alpha} \qquad (2.42) $$

This last equation is precisely the Kirchoff Current Law (KCL) which implies that $N-1$ cocycle fluxes are transient and vanish in the steady state, $\mathcal{J}_\gamma(t \to \infty) \to 0 \forall \gamma$. Therefore,
the analysis of the steady-state fluxes alone only requires the knowledge of a reduced number of degrees of freedom on the graph, whose number is given by the cardinality of cycles c.

**Forces** We now do the same for the affinity vector $\mathbf{A}(t) \in \mathcal{E}$ whose components are defined as:

$$A_e(t) = \log \left( \frac{r(e)p_{s(e)}(t)}{r(-e)p_{t(e)}(t)} \right) \qquad (2.43)$$

We put forward the following decomposition for the affinity vector:

$$\mathbf{A}(t) = \sum_{\gamma=1}^{N-1} \mathcal{A}_{\gamma}(t)\mathbf{c}^{\gamma} + \sum_{\alpha=N}^{E} \mathcal{A}_{\alpha}(t)\mathbf{e}^{\alpha} \qquad (2.44)$$

where coefficients $\mathcal{A}_{\gamma}(t) \equiv \mathbf{A}(t) \cdot \mathbf{e}^{\gamma}$ and $\mathcal{A}_{\alpha}(t) \equiv \mathbf{A}(t) \cdot \mathbf{c}^{\alpha}$, and we used the orthogonality properties Eqs. (2.13-2.12). Therefore, the coefficients $\mathcal{A}_{\alpha}$ for $N \le \alpha \le E$ corresponds to the cycle affinities, i.e. the affinities summed along cycles. Let us compute them explicitly by making use of Eq. (2.43):

$$\begin{aligned}
\forall \alpha, \mathcal{A}_{\alpha}(t) &= \mathbf{A}(t) \cdot \mathbf{c}^{\alpha} \\
&= \sum_{e} A_e(t) c_e^{\alpha} \\
&= \sum_{e} c_e^{\alpha} \log \left( \frac{r(e) p_{s(e)}(t)}{r(-e) p_{t(e)}(t)} \right) \qquad (2.45) \\
&= \log \left( \frac{\prod_{e} (r(e) p_{s(e)}(t))^{c_e^{\alpha}}}{\prod_{e} (r(-e) p_{t(e)}(t))^{c_e^{\alpha}}} \right) \\
&= \log \left( \frac{\prod_{e} (r(e))^{c_e^{\alpha}}}{\prod_{e} (r(-e))^{c_e^{\alpha}}} \right)
\end{aligned}$$

where the first equality was obtained using Eqs. (2.13) and (2.14). To go from the fourth to the fifth equality, we used the fact that the probabilities are vertex quantities and therefore their product along a cycle does not depend on the orientation. Several remarks:

* The cycle affinities only depend on the transition rates.

* If $\mathcal{A}_{\alpha} = 0 \forall \alpha$, then Eq. (2.45) ensures that the Kolmogorov criterion is fulfilled, namely that the product of transition rates across all cycles is the same in the forward and backward directions. This is a sufficient and necessary condition for the Markov chain to be reversible. It is also equivalent to the detailed balance property.

* The condition $\mathcal{A}_{\alpha} = 0 \forall \alpha$ is also called Kirchoff voltage law (KVL).

* For a reversible dynamics, i.e. when KVL is satisfied, all the affinities are conservative, which implies the existence of a potential vector $\mathbf{V}(t) \in \mathcal{V}$ such that:

$$\mathbf{A}(t) = \sum_{\gamma} \mathcal{A}_{\gamma}(t) \mathbf{c}^{\gamma} = -\mathbf{D}^{\top} \mathbf{V}(t) \qquad (2.46)$$

where we made use of the decomposition Eq. (2.44) and in the last equality we used the fact that cocycles span the $\mathrm{im}\,\mathbf{D}^{\top}$. Given that $\mathbf{D}^{\top}$ can be seen as a discrete gradient, Eq. (2.46) is a potential condition for conservative affinities.
* As a last remark, if fluxes fulfill KCL, namely $j(t) \in \ker D$, and affinities fulfill KVL, namely $A(t) \in \mathrm{im} D^\top$, then Tellegen's theorem is satisfied and we have:

$$j(t) \cdot A(t) = \sum_{\alpha, \gamma} \mathcal{J}_{\alpha} A_{\gamma} c^{\alpha} \cdot c^{\gamma} = 0 \qquad (2.47)$$

**Entropy production** Using the decomposition for the fluxes and for the affinities introduced above and the orthogonality relations expressed in subsection 2.2.2, we obtain a decomposition for the entropy production rate (EPR) as follows:

$$\begin{aligned}
\sigma(t) &\equiv A(t) \cdot j(t) \\
&= \left( \sum_{\gamma} A_{\gamma}(t) c^{\gamma} + \sum_{\alpha} A_{\alpha}(t) e^{\alpha} \right) \cdot \left( \sum_{\gamma'} \mathcal{J}_{\gamma'}(t) e^{\gamma'} + \sum_{\alpha'} \mathcal{J}_{\alpha'}(t) c^{\alpha'} \right) \\
&= \sum_{\gamma} \sum_{\gamma'} A_{\gamma}(t) \mathcal{J}_{\gamma'}(t) \underbrace{c^{\gamma} \cdot e^{\gamma'}}_{=\delta_{\gamma,\gamma'}} + \sum_{\gamma} \sum_{\alpha'} A_{\gamma}(t) \mathcal{J}_{\alpha'}(t) \underbrace{c^{\gamma} \cdot c^{\alpha'}}_{=0} \\
&\quad + \sum_{\alpha} \sum_{\gamma'} A_{\alpha}(t) \mathcal{J}_{\gamma'}(t) \underbrace{e^{\alpha} \cdot e^{\gamma'}}_{=0} + \sum_{\alpha} \sum_{\alpha'} A_{\alpha}(t) \mathcal{J}_{\alpha'}(t) \underbrace{e^{\alpha} \cdot c^{\alpha'}}_{\delta_{\alpha,\alpha'}} \\
&= \sum_{\gamma} \mathcal{J}_{\gamma}(t) A_{\gamma}(t) + \sum_{\alpha} \mathcal{J}_{\alpha}(t) A_{\alpha}(t)
\end{aligned} \qquad (2.48)$$

Thus there are again two cases:

* At steady-state, we know that the only remaining fluxes are the fluxes across the cycles $\mathcal{J}_{\alpha}(t)$. Therefore, in the long-time limit, the EPR simplifies to:

$$\sigma(\infty) = \sum_{\alpha} \mathcal{J}_{\alpha}(\infty) A_{\alpha}(\infty) \qquad (2.49)$$

* If the chain is reversible we have seen that $A_{\alpha} = 0 \forall \alpha$, therefore the EPR reduces to:

$$\sigma(t) = \sum_{\gamma} \mathcal{J}_{\gamma}(t) A_{\gamma}(t) \qquad (2.50)$$

and it vanishes in the long-time limit since $\mathcal{J}_{\gamma}(t \to \infty) \to 0 \forall \gamma$.

As a last remark, if the chain is both at steady state and reversible, i.e. if both KCL and KVL are fulfilled, then both terms in Eq. (2.48) vanish. The system is in an equilibrium steady-state with zero entropy production:

$$\sigma(t) = 0 \ \forall t \qquad (2.51)$$

## 2.5 An example of chemical reaction network

As a summary, we apply the Markov chain tree formula and the observable decomposition to the example of a chemical reaction network (CRN) inspired from biochemistry and shown in Figure 2.7. The reaction network describes the activation of a substrate $S \to S^*$ via the enzyme $E$ that binds to the substrate to form the complex $ES$. The latter complex also decomposes in the free enzyme and the new substrate $S^*$. A first observation is that the chemical network is non-linear and cannot be represented as a simple graph. Specifically, reactions 1 and 3 in panel a) of Figure 2.7 involve three species and cannot be represented as simple edges connecting a single vertex (species) to another one. This is
<figure><img src="image_33.png" /><figcaption>Figure 2.7: An example of CRN. a) Reaction 1 requires the enzyme E to bind to the substrate S to generate the complex ES. Similarly, reaction 3 requires the enzyme E to bind to the substrate S* to form the complex ES*. Consequently, the reaction network is not a graph and its dynamics is non linear. We need a hypergraph to represent it. b) By chemostating S and S*, the network becomes effectively linear. The concentration of S and S* are no longer dynamical variables and they enter in the dynamics via the rescaled transition rates $\tilde{k}_1^+ \equiv k_1^+ x_S^{chem}$ and $\tilde{k}_3^- x_{S^*}^{chem}$.</figcaption></figure>

typically the case in chemical reaction networks which involve interactions many-to-many and thus map to hypergraphs (we refer to [20] for a generalization of the graph-theoretical notions introduced here to hypergraphs). We assume in the following that the CRN follows the mass action kinetics and denote $x_i$ the concentration of species i. Using the mass action law, we can for instance write the kinetic equation for the concentration of ES through reaction 1 and 2. It reads:

$$ \frac{dx_{ES}}{dt} = k_1^+ x_S x_E - k_1^- x_{ES} - k_2^+ x_{ES} + k_2^- x_{ES^*} \quad (2.52) $$

Clearly, this kinetic equation is non linear due to the first term on the right-hand side. Here to illustrate the results developed in the previous subsections on graph theory, we choose to fix the concentration of S and S* via external chemostatting. The chemostating is performed by connecting the system with two reservoirs of respectively S and S*, held at concentrations $x_S^{chem}$ and $x_{S^*}^{chem}$. Then, the concentrations of S and S* no longer vary, and the network can be reduced effectively to a linear network on the states (E, ES, ES*), for which a graphical representation can be exploited (see Figure 2.7). Then, the kinetic equation (2.52) also reduces to a linear equation in the dynamical variables $x_E$ and $x_{ES}$:

$$ \frac{dx_{ES}}{dt} = \underbrace{k_1^+ x_S^{chem}}_{\tilde{k}_1^+} x_E - k_1^- x_{ES} - k_2^+ x_{ES} + k_2^- x_{ES^*} \quad (2.53) $$

We are thus left with a linear continuity equation formally analogous to the master equation (2.18) where the $x_i(t)$'s are the analogue of the probabilities $p_i(t)$ and the $k_i$'s are the analogue of the transition rates $r(i|j)$. Figure 2.7 panel b) represents the graph that can be obtained out of our linearized set of kinetic equations. For CRNs, the analogous of the
incidence matrix is called stoichiometric matrix. Each column of the stoichiometric matrix corresponds to a reaction of the CRN. We denote S the stoichiometric of our example which reads:

$$S = \begin{pmatrix} -1 & 0 & 1 \\ 1 & -1 & 0 \\ 0 & 1 & -1 \end{pmatrix} \qquad (2.54)$$

Since the graph of the CRN is unicyclic, the unique cycle vector reads:

$$c = \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix} \qquad (2.55)$$

where we dropped the index α. Since there is only one cycle, the force driving the CRN out of equilibrium is going to be the cycle affinity $A = c \cdot A$. We can compute it as follows:

$$\begin{aligned} A = A \cdot c &= \log \left( \frac{\tilde{k}_1^+ x_E k_2^+ x_{ES} k_3^+ x_{ES^*}}{k_1^- x_{ES} k_2 x_{ES^*} \tilde{k}_3^- x_E} \right) \\ &= \underbrace{\log \left( \frac{k_1^+ k_2^+ k_3^+}{k_1^- k_2^- k_3^-} \right)}_{1} + \underbrace{\log \left( \frac{x_S^{\text{chem}}}{x_{S^*}^{\text{chem}}} \right)}_{2} \end{aligned} \qquad (2.56)$$

We find two terms:

1. the first term depends only on the transition rates.

2. the second term depends only on the concentration of the reservoirs that we used to chemostat S and $S^*$.

The first term of equation (2.56) vanishes whenever the transition rates fulfill the Kolmogorov criterion which states that the following ratio must be one:

$$\frac{k_1^+ k_2^+ k_3^+}{k_1^- k_2^- k_3^-} = 1 \qquad (2.57)$$

The, A reduces to:

$$A = \log \left( \frac{x_S^{\text{chem}}}{x_{S^*}^{\text{chem}}} \right) \qquad (2.58)$$

and the out of equilibrium driving of our CRN relies only on the difference between the two chemostats connected to the CRN.

Finally, we can make the link with the Markov chain tree formula. Assuming Eq. (2.57) ensures that the chain is reversible and, in the absence of chemostatting, it relaxes to an equilibrium steady state. We then use a parameterization of the transition rates in terms of the standard chemical potential of thermodynamics $\mu_i^0$ of species i, namely:

$$\frac{k_e^+}{k_e^-} = e^{\mu_{s(e)}^0 - \mu_{t(e)}^0} \qquad (2.59)$$

This also ensures the following property at the level of spanning trees:

$$\frac{w(\mathcal{T}^i)}{w(\mathcal{T}^j)} = e^{\mu_j^0 - \mu_i^0} \qquad (2.60)$$
where $w(\mathcal{T}^i)$ is the product of the transition rates along the rooted spanning tree $\mathcal{T}^i$. Then, the Markov chain tree formula simplifies considerably leading to:

$$x_i(\infty) = \frac{w_i}{\sum_j w_j} = \frac{\sum_{\mathcal{T}} w(\mathcal{T}^i)}{\sum_j \sum_{\mathcal{T}} w(\mathcal{T}^j)} \qquad (2.61)$$

$$= \frac{\sum_{\mathcal{T}} w(\mathcal{T}^i)}{\sum_{\mathcal{T}} w(\mathcal{T}^i) \sum_j e^{\mu_i^0 - \mu_j^0}} = \frac{e^{-\mu_i^0}}{\sum_j e^{-\mu_j^0}} \qquad (2.62)$$

which is the equilibrium solution of thermodynamics.
# 3 Stochastic Thermodynamics

Massimiliano Esposito and Danilo Forastiere.<fnref n="6" /> These lecture notes contain an informal introduction to stochastic thermodynamics, which is a dynamical theory of mesoscopic systems in contact with external reservoirs which is fully compatible with statistical mechanics, and reproduces its main results when the full system is allowed to reach thermodynamic equilibrium. The stochastic dynamics of discrete systems is introduced in Section 3.1. For simplicity, in Section 3.2 we develop the stochastic dynamics and thermodynamics of a small system in contact with a single ideal reservoir in equilibrium, and at the average level. Then, in Section 3.3 we explain how to treat genuine nonequilibrium situations, in which the external reservoirs are characterized by different temperatures and therefore preclude the system from reaching equilibrium. Section 3.4 explains how to extend the formalism to cases in which the mesoscopic states of the system have internal entropy and how to deal with transfer of matter between reservoirs. Finally, Section 3.5 introduces the fundamental tools needed to develop stochastic thermodynamics at the level of the single stochastic trajectories, namely the fluctuation theorems.

## 3.1 Stochastic dynamics

We start from a system described by the Master Equation (with the graphical representation given in Fig. 3.1) with a possibly time-dependent generator W

$$d_t p_i = \sum_j W_{ij} p_j = \sum_{j \neq i} \underbrace{(W_{ij}(t)p_j - W_{ji}(t)p_i)}_{J_{ij}} \quad (3.1)$$

where conservation of normalization implies $\sum_j W_{ji} = 0$ or equivalently $W_{ii} = -\sum_{j \neq i} W_{ji}$, and where $J_{ij}$ denotes the probability current from j to i.

The steady state is obtained as right-null eigenvector of W and thus satisfies the linear equation $Wp^{ss} = 0$. It is unique for an ergodic Markov process<fnref n="7" /> by the Perron-Frobenius theorem [2,21]. If a time-dependent driving protocol is applied from the environment, the generator becomes time-dependent; in this case, we can define the instantaneous steady-state $p^{ss}(t)$ which satisfies $\sum_i W_{ji}(t)p_i^{ss}(t) = 0$ at every time and for each j, and which becomes time-independent when the time-dependent driving is stopped. In absence of time-dependent driving every initial condition relaxes to $p^{ss}$, as it can be seen introducing the relative entropy (or Kullback-Leibler divergence, [13,22]), which is the positive function defined as

$$D(p|p') \equiv \sum_i p_i \ln \frac{p_i}{p'_i} \quad (3.2)$$

$$= -\sum_i p_i \ln \frac{p'_i}{p_i} \geq -\sum_i p_i \left(1 - \frac{p'_i}{p_i}\right) = 0. \quad (3.3)$$

The inequality follows from the normalization $\sum_i p_i = \sum_i p'_i = 1$ and from the fact that $-\ln x \geq 1-x$ (which is easily obtained, for example, considering $\ln x = \int_1^x ds \frac{1}{s} \leq \int_1^x ds$ for $x > 1$ and analogously for $x < 1$, and which has the geometrical interpretation depicted in Fig. 3.2).

<footnote n="6">⁶ME was the lecturer; DF was the angel and wrote this chapter.</footnote>

<footnote n="7">⁷Defined by an irreducible generator W.</footnote>
<figure><img src="image_34.png" /><figcaption>Figure 3.1: Graphical representation of the ME (3.1). In a small time interval $dt$, a system in the state $i$ can jump to $j$ with a probability given by $W_{ji}dt$.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Figure 3.2: Graphical representation of the inequality $\ln x \le x - 1$.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>Figure 3.3: Closed system (no flux of matter) evolving at the fixed temperature T of the environment. The transition rates obey local detail balance in the form (3.7) since the system needs relax to equilibrium in absence of driving.</figcaption></figure>

The relative entropy $D(p|p^{ss})$ is zero only at the steady-state and its time-derivative is always negative, meaning that it is a Lyapunov function [23]. In fact

$$ -d_t D(p|p^{ss}(t)) = -\sum_i d_t p_i \ln \frac{p_i}{p_i^{ss}} - \underbrace{\sum_i p_i d_t \ln p_i}_{=0} + \sum_i p_i d_t \ln p_i^{ss}. \quad (3.4) $$

The second term vanishes by normalization using the ME (3.1). The third term only occurs in presence of time-dependent driving protocols. On the other hand, the first term obeys

$$ -\sum_i d_t p_i \ln \frac{p_i}{p_i^{ss}} = -\sum_{ij} W_{ij} p_j \ln \frac{p_i}{p_i^{ss}} = -\sum_{ij} W_{ij} p_j \ln \frac{p_i p_j^{ss}}{p_i^{ss} p_j} \quad (3.5) $$

$$ \geq \sum_{ij} W_{ij} p_j \left( 1 - \frac{p_i p_j^{ss}}{p_i^{ss} p_j} \right) = - \sum_i \frac{p_i}{p_i^{ss}} \sum_j W_{ij} p_j^{ss} = 0 \quad (3.6) $$

where we used the ME (3.1) in the first equality, the normalization property of the generator $\sum_i W_{ij} = 0$ in the second and fourth equalities, the inequality $-\ln \ge 1-x$ and finally the definition of the steady state. Therefore, when no driving is present and $d_t p^{ss} = 0$, we have the Lyapunov property $d_t D(p|p^{ss}) \le 0$.

As a side remark, notice that the matrix $\tilde{W}_{ij} \equiv W_{ij} \frac{p_i^{ss}}{p_i^{ss}}$ is a rate matrix ($\sum_i \tilde{W}_{ij} = 0$) which defines a different Master Equation rate matrix with the same steady state $p^{ss}$ as the original matrix $W$.

## 3.2 Stochastic thermodynamics for a single reservoir of energy

**Basic theory** We assign the energies of the states $\epsilon_i(t)$, the system is closed and in contact with a reservoir at temperature T (therefore $\beta = \frac{1}{k_B T}$), and we assume local detailed balance [24–26]:

$$ \frac{W_{ij}}{W_{ji}} = e^{-\beta(\epsilon_i - \epsilon_j)}. \quad (3.7) $$

Examples of rates satisfying eq. (3.7) are given later in this section, and are graphically depicted in Fig. 3.4.

Consider the average energy $\langle E \rangle = \sum_i \epsilon_i(t) p_i(t)$. Its rate of change in times is given by

$$ d_t \langle E \rangle = \underbrace{\sum_i p_i d_t \epsilon_i}_{\tilde{W}} + \underbrace{\sum_i \epsilon_i d_t p_i}_{\tilde{Q}} \quad (3.8) $$
This corresponds to the First Law of thermodynamics. The quantities $\dot{W}$ and $\dot{Q}$ have the dimensions of energy over time but they are not derivatives with respect to time for a generic protocol. A remark on notation: in these notes we use $d_t$, $\partial_t$ for total and partial derivatives and the overdot only to mean that a quantity has dimension of a rate of change, but without implying that it is a derivative.

The identification of the heat is motivated by the following identity

$$ \dot{Q} = \sum_{ij} W_{ij} p_j \epsilon_i = \sum_{ij} W_{ij} p_j (\epsilon_i - \epsilon_j) = \frac{1}{2} \sum_{ij} J_{ij} (\epsilon_i - \epsilon_j), \qquad (3.9) $$

obtained using the conservation of normalization $\sum_i W_{ij} = 0$ in the first equality and a symmetrization over nodes in the second. Inserting local detailed balance (3.7) in eq. (10.31) we get

$$ -\frac{\dot{Q}}{T} = \frac{k_B}{2} \sum_{ij} J_{ij} \ln \frac{W_{ij}}{W_{ji}}, \qquad (3.10) $$

from which we see that in a transition governed by LDB the entropy change in the reservoir is given by the heat divided by the temperature. The entropy of the system is

$$ S = -k_B \sum_i p_i \ln p_i. \qquad (3.11) $$

Its time derivative is the entropy change in the system and reads

$$ \begin{aligned} d_t S &= -k_B \sum_i d_t p_i \ln p_i = -k_B \sum_{ij} W_{ij} p_j (\ln p_i - \ln p_j) \\ &= -\frac{k_B}{2} \sum_{ij} J_{ij} \ln \frac{p_i}{p_j}. \end{aligned} \qquad (3.12) $$

Therefore the rate of entropy change in the universe (i.e. the entropy production rate) is obtained summing (3.10) and (3.12),

$$ \dot{\sigma} \equiv d_t S - \frac{\dot{Q}}{T} = \frac{k_B}{2} \sum_{ij} (W_{ij} p_j - W_{ji} p_i) \ln \frac{W_{ij} p_j}{W_{ji} p_i} \ge 0, \qquad (3.13) $$

where the inequality follows from $(a-b) \ln \frac{a}{b} \ge 0$ for positive $a, b$. The entropy production is zero only when $p$ satisfies detailed balance $W_{ij}p_j^{eq} = W_{ji}p_i^{eq}$. The equilibrium state is a steady state as $\sum_j W_{ij}p_j^{eq} = 0$ for which the canonical distribution holds. This can be obtained from the local detailed balance condition⁸ (3.7) by summing over $j$ and imposing normalization

$$ 1 = \sum_j p_j^{\text{eq}} = \sum_j \frac{W_{ji}}{W_{ij}} p_i^{\text{eq}} = p_i^{\text{eq}} e^{\beta \epsilon_i} \sum_j e^{-\beta \epsilon_j}, \qquad (3.14) $$

which implies

$$ p_i^{\text{eq}} = \frac{e^{-\beta \epsilon_i}}{\sum_j e^{-\beta \epsilon_j}} \equiv e^{-\beta(\epsilon_i - F^{\text{eq}})}, \qquad (3.15) $$

where we introduced the equilibrium free energy $F^{\text{eq}} \equiv -k_B T \ln \sum_j p_j^{\text{eq}}$. Since the system is in contact with a single reservoir, in absence of driving over energies it will relax to equilibrium.
<figure><img src="image_37.png" /><figcaption>Figure 3.4: Illustration of systems characterized by different types of rates. (a) Spin-boson (b) Arrhenius (c) Bose (d) Fermi (e) Mass-action chemistry</figcaption></figure>

**Examples of rates.** We list here some important examples of rates satisfying LDB. They characterize the systems depicted in Fig. 3.4.

**Spin boson.** The rates at which a spin in contact with a bosonic bath switches between up and down are given by

$$W_{+} = \frac{\hbar}{4|\Delta|} \gamma(|\Delta|) (\coth \frac{\beta\Delta}{2} - 1), \qquad (3.16a)$$

$$W_{-} = \frac{\hbar}{4|\Delta|} \gamma(|\Delta|) (\coth \frac{\beta\Delta}{2} + 1). \qquad (3.16b)$$

The LDB condition

$$\ln \frac{W^{+}}{W^{-}} = -\beta\Delta \qquad (3.17)$$

follows from the property $\coth \frac{x}{2} - 1 = \frac{2}{e^{x-1}} = e^{x}(\coth \frac{x}{2} + 1)$. Drude-Ullersma model [2,27] is obtained setting $\gamma(\Delta) = \frac{2}{\pi} \frac{\alpha^2}{\alpha^2 + \Delta^2}$.

**Arrhenius rates.** These rates can be obtained as the coarse-graining of a diffusion process in a multi-well potential. See e.g. [2] (Ch. XIII.6) and [26].

$$W_{ij} = \Gamma \exp \{-\beta(E_{ij}^* - \phi_j)\}, \qquad (3.18a)$$

$$W_{ji} = \Gamma \exp \{-\beta(E_{ji}^* - \phi_i)\}. \qquad (3.18b)$$

The rates involve the free energy $\phi_i = \epsilon_i - Ts_i$, where the entropic contribution arises from the curvature near the minima of the energy [see Fig. 3.4 (b)]. The barrier between

<footnote n="8">⁸For the present case of a system in equilibrium with a single reservoir, local detailed balance and detailed balance are equivalent (see e.g. [2] Ch.4).</footnote>
the two states is symmetric, $E_{ij}^* = E_{ji}^* = E^*$. The LDB (3.7) therefore takes the form

$$ \ln \frac{W_{ij}}{W_{ji}} = -\beta(\phi_i - \phi_j). \qquad (3.19) $$

**Fermi rates.** These rates describe the interaction of a two-level system with a bath obeying Fermi statistics. See [28].

Define $x = \epsilon - \mu$. The Fermi statistics is $f(x) = \frac{e^{-\beta x}}{e^{-\beta x+1}}$, therefore $1 - f(x) = \frac{1}{1+e^{-\beta x}}$.

The Fermi rates are given by

$$ W_+ = \Gamma f(\epsilon - \mu), \qquad (3.20a) $$

$$ W_- = \Gamma (1 - f(\epsilon - \mu)), \qquad (3.20b) $$

satisfying (3.7).

**Bose rates.** The Bose statistics is given by

$$ n(x) = \frac{e^{-x}}{1 - e^{-x}}. \qquad (3.21) $$

Correspondingly $1 + n(x) = \frac{1}{1-e^{-x}}$. The Bose rates are given by

$$ W_+ = \Gamma n(\epsilon - \mu), \qquad (3.22a) $$

$$ W_- = \Gamma (1 + n(\epsilon - \mu)), \qquad (3.22b) $$

satisfying (3.7).

**Mass-action chemistry.** Consider the chemical reaction $A + X \rightleftharpoons Y + B$ occurring in a container of volume $V$, where the concentrations of A and B are externally maintained constant by a chemostatting mechanism [29]. The mass-action transition rates $W_{(N_X,N_Y),(N'_X,N'_Y)}$ for going from the state with $(N'_X, N'_Y)$ particles of species X and Y respectively to the one with $(N_X, N_Y)$ particles are given by

$$ W_{(N_X-1,N_Y+1),(N_X,N_Y)} = k_+[A]N_X/V, \qquad (3.23a) $$

$$ W_{(N_X,N_Y),(N_X-1,N_Y+1)} = k_-[B](N_Y+1)/V. \qquad (3.23b) $$

The kinetic constants obey the constraint [29]

$$ \ln \frac{k_+}{k_-} = -\beta(\mu_A^0 + \mu_X^0 - \mu_B^0 - \mu_Y^0) \qquad (3.24) $$

where $\mu_i^0$ is the standard chemical potential of the species $i$ and for ideal solutions the chemical potential is obtained from the Gibbs free energy as $\mu_i = \frac{\partial}{\partial N_i} g_i(N_i) = \mu_i^0 + k_B T \ln N_i!$. It is then possible to show that (3.23) satisfy

$$ \ln \frac{W_{(N_X-1,N_Y+1),(N_X,N_Y)}}{W_{(N_X,N_Y),(N_X-1,N_Y+1)}} = -\beta\{g_X(N_X-1) + g(N_Y+1) - g(N_X) - g(N_Y)\}. \qquad (3.25) $$

Similar considerations can also be applied to non-ideal solutions [30].
**Nonequilibrium free energy** We can introduce the nonequilibrium free energy:

$$F = E - TS \qquad (3.26)$$

When the temperature of the environment is not allowed to vary, the nonequilibrium free energy changes in time according to

$$d_t F = d_t E - T d_t S, \qquad (3.27)$$

$$= \dot{W} - T\dot{\sigma}, \qquad (3.28)$$

using first (eq.(3.8)) and second law (eq. (3.13)) in the second equality.

Therefore we can express the entropy production rate in the form

$$\dot{\sigma} = \frac{\dot{W} - d_t F}{T}, \qquad (3.29)$$

which is the analogous of Kelvin's classic formulation of the second law of thermodynamics.

The nonequilibrium free energy can be related to the relative entropy between the actual state and the equilibrium state at a given temperature. In fact, using the definitions (3.2) and (3.26) we find

$$k_B T D(p|p^{\text{eq}}) = k_B T \sum_i p_i \ln p_i - k_B T \sum_i p_i \ln e^{-\beta(\epsilon_i - F^{\text{eq}})}, \qquad (3.30)$$

$$= -TS + E - F^{\text{eq}} = F - F^{\text{eq}}. \qquad (3.31)$$

It is important to note that in absence of driving (i.e. when no external power is provided, $\dot{W} = 0$) $F$ is a Lyapunov function for the system, meaning that dynamics will bring it to its minimum value (the equilibrium one). This can be seen using (3.29) combined with (3.31) gives

$$T\dot{\sigma} = \dot{W} - d_t F = -k_B T d_t D(p|p^{\text{eq}}) \geq 0. \qquad (3.32)$$

In the following, we will also make use of the identity

$$\dot{W} - d_t F^{\text{eq}} = -k_B T \sum_i p_i d_t \ln e^{-\beta(\epsilon_i - F^{\text{eq}})} = -k_B T \sum_i p_i d_t \ln p^{\text{eq}}. \qquad (3.33)$$

After summing and subtracting $d_t F^{\text{eq}}$ to the entropy production in eq. (3.29), eq. (3.33) gives

$$\dot{\sigma} = -k_B \left( \sum_i p_i d_t p_i^{\text{eq}} + d_t D(p|p^{\text{eq}}) \right) \qquad (3.34)$$

$$= -k_B \sum_i d_t p_i \ln \frac{p_i}{p_i^{\text{eq}}} = -k_B \sum_{ij} W_{ij} p_j \ln \frac{p_i}{p_i^{\text{eq}}}. \qquad (3.35)$$

**Special transformations**

**1. Reversible transformations.** If the transformation is extremely slow, every state of the system is at every time only slightly distant from its equilibrium value (see [31] for more details), $p_i = p_i^{\text{eq}} + \delta p_i$. Therefore

$$\ln \frac{p_i}{p_i^{\text{eq}}} = \ln \left( 1 + \frac{\delta p_i}{p_i^{\text{eq}}} \right) \approx \frac{\delta p_i}{p_i^{\text{eq}}} - \frac{1}{2} \left( \frac{\delta p_i}{p_i^{\text{eq}}} \right)^2 + O \left( \left( \frac{\delta p_i}{p_i^{\text{eq}}} \right)^3 \right). \qquad (3.36)$$
The resulting entropy production is a second order quantity in the relative variation of the states,

$$ \dot{\sigma} \approx -k_B \sum_{ij} W_{ij} p_j^{\text{eq}} \left( 1 + \frac{\delta p_j}{p_j^{\text{eq}}} \right) \left( \frac{\delta p_i}{p_i^{\text{eq}}} + \frac{1}{2} \left( \frac{\delta p_i}{p_i^{\text{eq}}} \right)^2 \right) \quad (3.37) $$

$$ = -k_B \sum_{ij} W_{ij} p_j^{\text{eq}} \frac{\delta p_i}{p_i^{\text{eq}}} + k_B \sum_{ij} W_{ij} p_j^{\text{eq}} \frac{\delta p_i}{p_i^{\text{eq}}} \frac{\delta p_j}{p_j^{\text{eq}}} - \frac{k_B}{2} \sum_{ij} W_{ij} p_j^{\text{eq}} \left( \frac{\delta p_i}{p_i^{\text{eq}}} \right)^2, \quad (3.38) $$

since the first sum vanishes by normalization.

**2. Sudden switch.** Driving faster than the typical relaxation rate of the dynamics does not give time to p to change in time. Therefore, like in an isolated system

$$ \dot{\sigma} \approx 0, \dot{Q} \approx 0, d_t S \approx 0. \quad (3.39) $$

Crucially, however, $\dot{W} = d_t F \neq 0$.

**Nonequilibrium state as a resource** We consider again the dissipated free energy expressed via eq. (3.29)

$$ T \dot{\sigma} = \dot{W} - d_t F^{\text{eq}} - k_B T d_t D(p|p^{\text{eq}}) \geq 0 \quad (3.40) $$

Integrating it with respect to time and defining the total entropy production $\sigma \equiv \int_0^t dt' \dot{\sigma}(t')$ we obtain the integrated form of the Kelvin formulation of the second law

$$ T\sigma = \underbrace{W - \Delta F^{\text{eq}}}_{W_{\text{irr}}} - k_B T D(p(t)|p^{\text{eq}}(t)) + k_B T D(p(0)|p^{\text{eq}}(0)) \geq 0, \quad (3.41) $$

where the equilibrium distribution at different time can change because the energy levels $\{\epsilon_i\}$ are varied by the externally supplied work. Note, however, that the above calculation is restricted to isothermal conditions.

In a transformation between equilibrium states ($p(0) = p^{\text{eq}}(0)$ and $p(t) = p^{\text{eq}}(t)$), the contribution coming from the relative entropy cancels and eq. (3.41) implies that irreversible work is positive ($W_{\text{irr}} \geq 0$), or equivalently

$$ W \geq \Delta F^{\text{eq}}. \quad (3.42) $$

This is exactly Kelvin's formulation of the second law for transformations between equilibrium states.

However, preparing the system in a nonequilibrium state $p(0) \neq p^{\text{eq}}(0)$ allows the irreversible work

$$ W_{\text{irr}} = \underbrace{T\sigma}_{\geq 0} + \underbrace{k_B T D(p(t)|p^{\text{eq}}(t))}_{\geq 0} - \underbrace{k_B T D(p(0)|p^{\text{eq}}(0))}_{\leq 0} \quad (3.43) $$

to become negative due to the nonequilibrium free energy stored in the initial condition, thus achieving extraction from the nonequilibrium initial state, which acts as a resource for the environment [32] (see Fig. 3.5).
<figure><img src="image_38.png" /><figcaption>Figure 3.5: Extreme cases for free energy transduction, with red crosses corresponding to vanishing contributions. The resource term $k_B T D(P_0 | P_0^{\text{eq}})$ due to the nonequilibrium initial state can be either completely lost ($W = 0$ and $\sigma = k_B D(P_0 | P_0^{\text{eq}})$) if no extraction protocol is applied, or optimally extracted via a reversible, quasi-static protocol (for which $\sigma = 0$ and $W = -k_B T D(P_0 | P_0^{\text{eq}})$) following a sudden change of the potential from two to single well, followed by an adiabatic deformation back to two wells.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Figure 3.6: Systems in contact with multiple reservoirs. (a) Schematic depiction of a generic system (b) scheme of a 2-levels system in contact with two baths at different temperatures.</figcaption></figure>
## 3.3 ST for multiple reservoirs

We can generalize the above treatment to the important case of multiple reservoirs.

The local detailed balance takes the form

$$ \ln \frac{W_{ij}^{(\nu)}}{W_{ji}^{(\nu)}} = -\beta_{\nu}(\epsilon_i - \epsilon_j), \qquad (3.44) $$

where the index $\nu$ now spans the different reservoirs in contact with the system. Each of the reservoirs, in isolation, would impose a well defined equilibrium state as the steady state of the dynamics of the system. However, the simultaneous interaction of the system with different reservoirs frustrates the relaxation to any of the equilibria identified by each of them, and consequently heat and particle flows ensue. The dynamics follows the master equation (3.1), where we assume that the reservoirs affect the system independently of each other, leading to an additive generator $W = \sum_{\nu} W^{(\nu)}$ [24].

We identify the heat in each reservoir as being

$$ \dot{Q}^{(\nu)} = \frac{1}{2} \sum_{ij} \frac{1}{2} (W_{ij}^{(\nu)} p_j - W_{ji}^{(\nu)} p_i) \frac{W_{ij}^{(\nu)}}{W_{ji}^{(\nu)}} \qquad (3.45) $$

This gives for the entropy production rate appearing in the second law (3.13)

$$ \dot{\sigma} = d_t S - \sum_{\nu} \frac{\dot{Q}^{(\nu)}}{T} \qquad (3.46) $$

$$ = k_B \sum_{\nu, i, j} \frac{1}{2} (W_{ij}^{(\nu)} p_j - W_{ji}^{(\nu)} p_i) \ln \frac{W_{ij}^{(\nu)} p_j}{W_{ji}^{(\nu)} p_i} \geq 0, \qquad (3.47) $$

where the inequality follows from $(a-b) \ln \frac{a}{b} \geq 0$ for positive $a, b$.

Applying the log-sum inequality [22] $\sum_n a_n \ln \frac{a_n}{b_n} \geq (\sum a_n) \ln \frac{\sum_n a_n}{\sum_n b_n}$ (for nonnegative $a_n, b_n$), to the sum over the index $\nu$ of the reservoirs in eq. (3.47), after defining the probability flow $j_{ij}^{(\nu)} = W_{ij}^{(\nu)} p_j$ due to each of the reservoirs, we have

$$ \dot{\sigma} = \frac{1}{2} \sum_{\nu, i, j} (j_{ij}^{\nu} - j_{ji}^{\nu}) \ln \frac{j_{ij}^{\nu}}{j_{ji}^{\nu}} = \sum_{\nu, i, j} j_{ij}^{\nu} \ln \frac{j_{ij}^{\nu}}{j_{ji}^{\nu}} \geq \sum_{i, j} \left( \sum_{\nu} W_{ij}^{\nu} \right) p_j \ln \frac{\left( \sum_{\nu} W_{ij}^{\nu} \right) p_i}{\left( \sum_{\nu} W_{ji}^{\nu} \right) p_j}. \qquad (3.48) $$

This means that lumping together rates coming from the interaction with different reservoirs underestimates the entropy production. In particular, there could be detailed balance effective models that do not resolve between different reservoirs and that therefore estimate a vanishing entropy production even for a system which is out of equilibrium.

**Example:** The above fact can be seen easily by considering the 2-levels system ($p_0 = 1 - p_1$) in contact with two reservoirs of Fig. 3.6 (b) which is governed by the master equation $d_t p_0 = W_{01}p_1 - W_{10}p_0$. It reaches a stationary state when the detailed balance condition $W_{01}p_1 = W_{10}p_0$ is satisfied, thus leading to an estimate of zero entropy production when the latter is computed with the lumped rates $W$. However, since $W = \sum_{\nu=l,r} W^{(\nu)}$, the steady state entropy production rate computed using the correct formula (3.47) is nonzero whenever there is a difference in the temperatures of the two reservoirs. Therefore, when in contact with multiple reservoirs the steady state toward which the state of the systems relaxes by virtue of the Perron-Frobenius theorem will not coincide with the equilibrium state given by the Gibbs distribution (3.15). Therefore, in
<figure><img src="image_40.png" /><figcaption>Figure 3.7: System coupled to energy and matter reservoirs.</figcaption></figure>

absence of external time-dependent driving, after a sufficiently long time the system will be in a nonequilibrium steady state (NESS) characterized by a positive entropy production rate $\dot{\sigma}^{ss} \ge 0$.

The first law now takes the form — obtained following the same steps as for eq. (10.31)

$$d_t E = \sum_i d_t \epsilon p_i + \sum_i \epsilon d_t p_i = \dot{W} + \sum_{\nu, i, j} j_{ij}^{\nu} (\epsilon_i - \epsilon_j) \quad (3.49)$$

$$= \dot{W} + \sum_{\nu, i, j} \frac{1}{2} (j_{ij}^{\nu} - j_{ji}^{\nu}) (\epsilon_i - \epsilon_j) = \dot{W} + \sum_{\nu} \dot{Q}^{\nu}, \quad (3.50)$$

where $\dot{Q}^{\nu}$ is the heat flowing in the $\nu$-th reservoir.

If we are in presence of reservoirs indexed by $\nu = 0, \dots, N$, and we choose $T_0$ as the reference temperature which appears in the nonequilibrium free energy definition (3.26), we can rewrite the entropy production rate (3.46)

$$T_0 \dot{\sigma} = T_0 d_t S - \dot{Q}^0 - \sum_{\nu \ge 1} \frac{T_0}{T_\nu} Q^\nu \quad (3.51)$$

$$= \dot{W} - d_t F^0 + \sum_{\nu \ge 1} \left(1 - \frac{T_0}{T_\nu}\right) Q^\nu \quad (3.52)$$

$$= \dot{W} - d_t F^{\text{eq},0} + \sum_{\nu \ge 1} \left(1 - \frac{T_0}{T_\nu}\right) Q^\nu - k_B T_0 d_t D(p(t)|p^{\text{eq}}(t)). \quad (3.53)$$

Two remarks are needed here. First, the second equality we used the first law (3.50) to solve for $\dot{Q}^0$, the definition of nonequilibrium free energy (3.26) and we have identified the Carnot efficiency $\eta_C^\nu = 1 - \frac{T_\nu}{T_0}$ of the heat exchanged with the $\nu$-th reservoir, which is zero when evaluated for $T_0$ (i.e., no power output can be extracted from a single heat reservoir, that is from an equilibrium environment). Second, in the third equality we made use of the identity for the nonequilibrium free energy in terms of the equilibrium one plus a relative entropy, eq. (3.31).

In a nonequilibrium steady state it reduces to

$$\dot{\sigma}^{ss} = \sum_{\nu} \left( \frac{1}{T_0} - \frac{1}{T_{\nu}} \right) \dot{Q}^{\nu}, \quad (3.54)$$

which is analogous to the expression obtained in macroscopic nonequilibrium thermodynamics [33].

## 3.4 Internal entropy and matter transfer

We can generalize the description to include states with internal entropies and to allow for the exchange of particles with the reservoirs.
Each node is now characterized by a number of particles $n_i$, a mesoscopic energy $\epsilon_i = e_i n_i$ ($e_i$ being the average energy per particle in state $i$), and an internal entropy $s_i$. The average values of these observables in a state characterized by the distribution $p$ are respectively $\langle E \rangle = \sum_i \epsilon_i p_i$, $\langle N \rangle = \sum_i n_i p_i$ and $S = \sum_i (s_i - \ln p_i)p_i$. Notice that the Shannon entropy has been slightly extended to take into account the mesoscopic entropy. Energy and particles are conserved quantities for the (isolated) total system obtained considering the mesoscopic subsystem together with all the reservoirs. Given a reservoir with temperature $T^\nu$ and chemical potential $\mu^\nu$, we can define the nonequilibrium free energy associated with the mesoscopic states as $\phi_i = \epsilon_i - \mu^\nu n_i - k_B T^\nu s_i$. The local detailed balance becomes

$$ \ln \frac{W_{ij}^{\nu}}{W_{ji}^{\nu}} = -\beta^{\nu} (\phi_{i}^{\nu} - \phi_{j}^{\nu}) . \qquad (3.55) $$

The fluxes of the conserved quantities can be derived considering their balance equations. For the number of particles we have (repeating the same steps used for eq. (10.31))

$$ d_t \langle N \rangle = \sum_{\nu ij} W_{ij}^{\nu} p_j n_i = \sum_{\nu ij} \frac{1}{2} (W_{ij}^{\nu} p_j - W_{ji}^{\nu} p_i) (n_i - n_j) = \sum_{\nu} I_N^{\nu}, \qquad (3.56) $$

introducing the particle currents for each reservoir, $I_N^\nu = \sum_{ij} \frac{1}{2}(W_{ij}^\nu p_j - W_{ji}^\nu p_i)(n_i - n_j)$. Analogously, we introduce the energy flows $I_E^\nu = \frac{1}{2} \sum_{ij} (W_{ij}^\nu p_j - W_{ji}^\nu p_i)(\epsilon_i - \epsilon_j)$ to write the first law in a form in which the chemical work contribution $\dot{W}_{\text{chem}}$ appears explicitly,

$$ d_t \langle E \rangle = \underbrace{\sum_i d_t \epsilon_i p_i}_{\dot{W}_d} + \underbrace{\sum_{\nu} (I_E^{\nu} - \mu^{\nu} I_N^{\nu})}_{\dot{Q}^{\nu}} + \underbrace{\sum_{\nu} \mu^{\nu} I_N^{\nu}}_{\dot{W}_{\text{chem}}}. \qquad (3.57) $$

The identification of heat flow with

$$ \dot{Q}^{\nu} = \frac{1}{2} \sum_{ij} (W_{ij}^{\nu} p_j - W_{ji}^{\nu} p_i) ((\epsilon_i - \mu^{\nu} n_i) - (\epsilon_j - \mu^{\nu} n_j)) \qquad (3.58) $$

$$ = \frac{1}{2} \sum_{ij} (W_{ij}^{\nu} p_j - W_{ji}^{\nu} p_i) ((\phi_i + T^{\nu} s_i) - (\phi_j + T^{\nu} s_j)) \qquad (3.59) $$

takes into account the fact that part of the energy flow is associated to the transfer of particles between different reservoirs, and therefore it should not be counted as heat (as it does not “heat” the reservoir). Indeed, one can easily verify that at equilibrium with a single reservoir eq. (3.59) gives the Clausius formula for the entropy change for a given heat flux once the equilibrium distribution $p_i^{\text{eq}} = Z^{-1}e^{-\beta\phi_i}$ is introduced.

The second law can be obtained computing

$$ d_t S = \sum_i d_t p_i (s_i - \ln p_i) = \sum_{\nu, i, j} \frac{1}{2} (W_{ij}^{\nu} p_i - W_{ji}^{\nu} p_i) \left(s_i - s_j - \ln \frac{p_i}{p_j}\right). \qquad (3.60) $$

The expression for the entropy production rate, encoding the second law, is consequently

$$ \dot{\sigma} = d_t S - \sum_{\nu} \frac{\dot{Q}^{\nu}}{T^{\nu}} = \sum_{\nu, i, j} \frac{1}{2} (W_{ij}^{\nu} p_j - W_{ji}^{\nu} p_i) \ln \frac{W_{ij}^{\nu} p_i}{W_{ji}^{\nu} p_i} \geq 0. \qquad (3.61) $$
Alternatively, using $\dot{Q}^{\nu} = I_{E}^{\nu} - \mu^{\nu} I_{N}^{\nu}$ and using the first law (3.57) to isolate the reservoir with $\nu = 0$, used as a reference, we can write

$$T_0 \dot{\sigma} = T_0 d_t S - I_E^0 + \mu^0 I_N^0 - \sum_{\nu \ge 1} (I_E^\nu - \mu^\nu I_N^\nu) \quad (3.62)$$

$$= \dot{W}_d - (d_t \langle E \rangle - \mu_0 d_t \langle N \rangle - T_0 d_t S) + \sum_{\nu} \left( \frac{1}{T_0} - \frac{1}{T_{\nu}} \right) I_E^{\nu} - \sum_{\nu} \left( \frac{\mu_{\nu}}{T_{\nu}} - \frac{\mu_0}{T_0} \right) I_N^{\nu} \quad (3.63)$$

and therefore

$$\dot{\sigma} = \frac{\dot{W}_d - d_t \Phi^0}{T_0} + \sum_{\nu \ge 1} (F_E^{\nu} I_E^{\nu} + F_N^{\nu} I_N^{\nu}) \quad (3.64)$$

with the nonconservative forces defined as $F_E^{\nu} = \frac{1}{T_0} - \frac{1}{T_{\nu}}$ and $F_N^{\nu} = \frac{\mu_{\nu}}{T_{\nu}} - \frac{\mu_0}{T_0}$. The decomposition of the entropy production (3.64) obtained here can be refined using the conceptually similar but more general framework of Refs. [12, 13]. There, it is also proved that the number of fundamental force $N_F$ can be expressed as the difference between the number of intensive variables associated to the reservoirs $N_I$ and the number of conserved quantities that characterize the total system $N_C$, namely

$$N_F = N_I - N_C. \quad (3.65)$$

The important case of tight-coupling between currents, namely when the transfer process for different conserved quantities (e.g. the one of matter and energy) occurs simultaneously, is discussed in Ref. [34].

## 3.5 Fluctuation theorems and thermodynamic uncertainty relations

Until now, we have considered the average value of thermodynamic observables, where the averaging process sums over all possible stochastic trajectories compatible with the dynamics. Stochastic thermodynamics, however, can be formulated also at the level of the single trajectories [35,36]. Let's indicate with $\Gamma$ a trajectory in the state space. We define an operation of time reversal $\Gamma \to \tilde{\Gamma}$, which changes the order of the visited states and the sign of all the observable which are odd under time-reversal. Furthermore, it reverses the time-dependence contained in the time-dependent driving protocols, if present.

We have the following trajectory level entropy production

$$\sigma(\Gamma) = k_B \ln \frac{P(\Gamma)}{\tilde{P}(\Gamma)} \quad (3.66)$$

which is a trajectory-wise version of the local detailed balance condition (3.7) and is related to the so-called detailed fluctuation theorem [37–41].

An important result is the integral fluctuation theorem [42], that in this context is an immediate consequence of the rearrangement of eq. (3.66)

$$\langle e^{-\frac{\sigma}{k_B}} \rangle = \sum_{\Gamma} P(\Gamma) \frac{\tilde{P}(\tilde{\Gamma})}{P(\Gamma)} = \sum_{\Gamma} \tilde{P}(\tilde{\Gamma}) = \sum_{\tilde{\Gamma}} \tilde{P}(\tilde{\Gamma}) = 1. \quad (3.67)$$

The third equality is obtained recognizing that the sum over all paths can be enumerated equivalently using the forward or the reversed trajectories.
The average of (3.66) gives the entropy production at the average level

$$ \langle \sigma \rangle = k_B \sum_{\Gamma} P(\Gamma) \ln \frac{P(\Gamma)}{\bar{P}(\Gamma)} \ge 0. \qquad (3.68) $$

The inequality is a consequence of Jensen's inequality, as $1 = \langle e^{-\sigma/k_B} \rangle \ge e^{\langle \sigma \rangle/k_B} \ge 1 - \frac{\langle \sigma \rangle}{k_B}$.

Another important class of results concerning fluctuations is the one of thermodynamic uncertainty relations (TUR). TURs are relations which bound the signal-to-noise ratio of a current $J$ (or another measure of the precision of a signal) in terms of the entropy production needed to achieve it. The original TUR [43] states that

$$ \mathrm{pr}(J) \equiv \frac{\langle J \rangle^2}{\mathrm{Var}(J)} \le \frac{\langle \sigma \rangle}{2k_B}, \qquad (3.69) $$

and therefore provides a bound on a dynamical characterization of a stochastic system (the precision) in terms of a thermodynamic observable (the EPR), which can be exploited either to optimize the former or to infer the latter. For a general perspective unifying different types of uncertainty relations, see [44].
# 4 Deterministic Chemical Reaction Networks

Francesco Avanzini, Shesha Gopal Marehalli Srinivas, Emanuele Penocchio and Massimiliano Esposito.<fnref n="9" /> We formulate a nonequilibrium thermodynamic theory for open chemical reaction networks (CRNs) described by deterministic rate equations following the law of mass-action. The conservation laws of CRNs are used to decompose the entropy production into a potential change and two work contributions. One work contribution is due to time-dependent changes in the externally controlled chemostat concentrations. The other accounts for the flows maintained across the CRN by nonconservative forces breaking the detailed balance condition.

## 4.1 Introduction

Open CRNs constitute the underlying mechanism of many complex processes in biosystems [45] and synthetic chemistry [46] across many scales. Prototypical examples are information processing [47], oscillations [48], self-replication [49–51], self-assembly [52, 53], and molecular machines [54, 55]. These processes operate out of equilibrium: nonzero reaction currents are sustained by continuously harvesting free energy from environment (represented in terms of reservoirs of chemical species called chemostats) [56]. Their energetics can be characterized on rigorous grounds using a nonequilibrium thermodynamic theory for CRNs. Such a theory has been developed in recent years for CRNs undergoing a stochastic [29, 57, 58] or a deterministic dynamics [10, 30, 59–61]. This theory has been used to quantify the energetic costs of maintaining coherent oscillations [62] and sustaining growth of copolymers [63–65] and biomolecules [66]; the efficiency of central metabolism [67] in prokaryotes; the internal information transfer in a model of chemically-driven self-assembly and an experimental light-driven bimolecular motor [68]. It was also formulated for reaction diffusion systems, where chemical reactions can be used to create patterns [69, 70] and waves [71].

In these lecture notes, we focus on deterministic CRNs whose dynamics follows the law of mass-action. In particular, we re-derive, in a self-contained way, the formulation of nonequilibrium thermodynamics for CRNs developed in Refs. [29, 30, 60]. The lecture notes are organized as follows. We introduce the basic setup for the description of CRNs in Sec. 4.2 and discuss their dynamics in Sec. 4.3. We then build nonequilibrium thermodynamics on top of the CRN dynamics. We start in Sec. 4.4 where we introduce the notion of thermodynamic consistency: closed CRNs must be detailed balanced, namely, they must relax towards an equilibrium state. This allows us to derive a condition known in stochastic thermodynamics as local detailed balance establishing a correspondence between dynamic (i.e., reaction fluxes) and thermodynamic (i.e., chemical potentials) quantities. In Sec. 4.5, we use the local detailed balance to obtain the nonequilibrium formulation of the first and second law of thermodynamics for closed CRNs. In Sec. 4.6, we derive the first and second law of thermodynamics for open CRNs. Crucially, we show how to use the conservation laws of CRNs (defined in Subs. 4.3.3) to rewrite the second law in such a way as to split the free energy exchanged with the chemostats into two work contributions. One contribution, named driving work, accounts for the free energy exchanged with the chemostats via a time dependent manipulation of CRNs. The other contribution, named nonconservative work, accounts for the nonconservative forces created via the exchanges with the chemostats that break the detailed balance condition and can maintain CRNs out of equilibrium. This decomposition of the second law is a major result as it identifies

<footnote n="9">⁹FA was the lecturer and wrote this chapter; SGMS was the angel; EP and ME contributed to the selection and organization of the chapter content.</footnote>
the specific mechanism, in terms of nonconservative forces, that can maintain CRNs out of equilibrium.

**Disclaimer.** These lecture notes represent a first (partial) summary of the recent developments in nonequilibrium thermodynamics of CRNs that Francesco Avanzini, Mas-similiano Esposito and Emanuele Penocchio (authors are listed in alphabetic order) are planning to review more extensively in another contribution. The course “Thermodynamics of Deterministic Chemical Reaction Networks” of the school (Post)Modern Thermodynamics (5-9 December 2022, Luxembourg) was given by Francesco Avanzini and prepared together with Shesha Gopal Marehalli Srinivas.

## 4.2 Basic Setup

A CRN is defined as a set of chemical species {$Z_\alpha$}, identified by $\alpha \in S = \{1, 2, \dots, n_s\}$, that are interconverted via chemical reactions, identified by the index $\rho \in \mathcal{R} = \{1, 2, \dots, n_r\}$. Each reaction $\rho$ is assumed here to be reversible and represented by an equation like

$$ \sum_{\alpha \in S} \nu_{\alpha, \rho}^{+} Z_{\alpha} \underset{\rho}{\stackrel{\rho}{\rightleftharpoons}} \sum_{\alpha \in S} \nu_{\alpha, \rho}^{-} Z_{\alpha}, \qquad (4.1) $$

where $\nu_{\alpha, \rho}^{+}$ (resp. $\nu_{\alpha, \rho}^{-}$) is the stoichiometric coefficients of $Z_\alpha$ in the forward (resp. backward) transformation. The net stoichiometry is encoded in the $(n_s \times n_r)$ stoichiometric matrix $S$ whose columns are defined as

$$ S_{\rho} := \nu_{\rho}^{-} - \nu_{\rho}^{+}, \qquad (4.2) $$

with $\nu_{\rho}^{\pm} := (\dots, \nu_{\alpha, \rho}^{\pm}, \dots)^{\mathrm{T}}_{\alpha \in S}$.

*Physical Remark.* We consider here systems, known as ideal dilute solutions, where the $n_s$ species {$Z_\alpha$} are mixed together with a non-reacting, very abundant species $Z_0$ called solvent. The solvent maintains the temperature $T$ and the volume $V$ of the system constant.

*Example.* We now consider the following CRN,

$$ \begin{array}{l} F + E \underset{EF}{\overset{1}{\rightleftharpoons}} EF \\ \\ EF \underset{EW}{\overset{2}{\rightleftharpoons}} EW \\ \\ EW \underset{E+W}{\overset{3}{\rightleftharpoons}} E + W \qquad (4.3) \\ \\ S + EF \underset{EFS}{\overset{4}{\rightleftharpoons}} EFS \\ \\ EFS \underset{EW+P}{\overset{5}{\rightleftharpoons}} EW + P \end{array} $$

representing a minimal metabolic process. Indeed, the interconversion of the substrate S (representing for instance ADP) into the product P (representing for instance ATP) is catalyzed by the enzyme E and powered by the interconversion of the fuel F (representing for instance nutrients) into the waste W (representing for instance $CO_2$). The species EF and EW are the complexes enzyme-fuel and enzyme-waste, respectively. The species EFS is the complex obtained by binding both the fuel and the substrate to the enzyme.
The stoichiometric matrix of the CRN (4.3) reads

$$ S = \begin{array}{c} \begin{array}{c} \begin{array}{ccccc} 1 & 2 & 3 & 4 & 5 \\ E & EF & EW & EFS & P \\ W & S & F \end{array} \end{array} \begin{array}{ccccc} -1 & 0 & 1 & 0 & 0 \\ 1 & -1 & 0 & -1 & 0 \\ 0 & 1 & -1 & 0 & 1 \\ 0 & 0 & 0 & 1 & -1 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & -1 & 0 \\ -1 & 0 & 0 & 0 & 0 \end{array} \end{array}, \quad (4.4) $$

where, for instance, the first column specifies that 1 molecule of E and 1 molecule of F are consumed every time the first reaction occurs (in the forward direction), while 1 molecule of EF is produced.

*Mathematical Remark.* The same stoichiometric matrix might correspond to different CRNs. Consider for instance two CRNs composed of the three species S, P, and E. In one CRN, the species undergo the chemical reaction

$$ S \rightleftharpoons P. \quad (4.5) $$

In the other CRN, the species undergo the chemical reaction

$$ S + E \rightleftharpoons P + E. \quad (4.6) $$

According to Eq. (4.2), both CRNs have the same stoichiometric matrix.

**Disclaimer.** The CRN (4.3) has been chosen to illustrate the theory discussed in these lecture notes and does not necessarily represent any realistic chemical process.

## 4.3 Dynamics

The evolution of deterministic CRNs with constant volume is specified by the (vector of the) concentrations of chemical species: $[Z] = (..., [Z_\alpha], ... )_{\alpha \in S}$.

*Notation Remark.* Throughout these notes, we will omit for compactness of notation the time dependence of any quantity (e.g., concentrations and thermodynamic fluxes, like heat, entropy flow, and entropy production rate).

### 4.3.1 Rate Equation

In open CRNs, the concentration vector follows the rate equation

$$ \frac{d}{dt}[Z] = S\mathbf{J}([Z]) + \mathbf{I}. \quad (4.7) $$

where $\mathbf{J}([Z])$ is the reaction current vector and $\mathbf{I}$ is the exchange current vector.

Each entry $J_\rho$ of $\mathbf{J}$ quantifies the net current of reaction $\rho$. If $J_\rho > 0$ (resp. $J_\rho < 0$), reaction $\rho$ occurs in the forward (resp. backward) direction, namely, from left to right (resp. from right to left) in Eq. (4.1). Each net current is given by the difference between the forward and backward fluxes,

$$ J_{\rho}([Z]) = r_{\rho}^{+}([Z]) - r_{\rho}^{-}([Z]), \quad (4.8) $$
which satisfy the law of mass-action, i.e., the fluxes are proportional to the concentrations of the reactants to the power of their stoichiometric coefficients:

$$r_{\rho}^{\pm}([Z]) = k_{\rho}^{\pm}[Z]^{\nu_{\rho}^{\pm}}, \qquad (4.9)$$

with $[Z]^{\nu_{\rho}^{\pm}} = \prod_{\alpha \in S}[Z_{\alpha}]^{\nu_{\alpha, \rho}^{\pm}}$ and $k_{\rho}^{\pm}$ being the so-called kinetic constants. This is physically justified for ideal dilute solutions: interactions between chemical species are negligible and the solvent is much more abundant, i.e., $[Z_0] \gg \sum_{\alpha \in S}[Z_{\alpha}]$. Therefore, reacting species behave as an ideal gas mixture reacting at rates that are proportional to their concentrations.

Each entry $I_{\alpha}$ of $I$ quantifies the net current at which the species $Z_{\alpha}$ is exchanged with the environment. This can represent for instance i) exchanges with external particle reservoirs known as chemostats, or ii) the net effect of other non-specified chemical reactions.

We can now split the set of chemical species $\{Z_{\alpha}\}$, and of the corresponding set of indexes $S$, into two disjoint sets by using the rate equation (4.7). The $n_X$ species $\{X_{\alpha}\}$ (with $\alpha \in S_X$) are not exchanged with the environment, i.e., $I_{\alpha} = 0$ for the whole dynamics, and are thus called *internal*. The $n_Y$ species $\{Y_{\alpha}\}$ (with $\alpha \in S_Y$) are exchanged with the environment, i.e., $I_{\alpha} \neq 0$ for some moments of the dynamics, and are thus called *exchanged*. By applying the same splitting to the concentration vector

$$[Z] = ([X], [Y])^{\mathrm{T}}, \qquad (4.10)$$

and the stoichiometric matrix

$$S = \begin{pmatrix} S^X \\ S^Y \end{pmatrix}, \qquad (4.11)$$

the rate equation (4.7) can be rewritten as

$$\frac{d}{dt}[X] = S^X J([X], [Y]), \qquad (4.12)$$

$$\frac{d}{dt}[Y] = S^Y J([X], [Y]) + I^Y, \qquad (4.13)$$

where $I^Y$ collects the non-null entries of $I$.

*Remark.* In these lecture notes, we consider CRNs coupled with chemostats (the $S_Y$ species are thus said to be chemostatted). This means that the concentrations $[Y]$ are not dynamical variables, but are controlled by the external chemostats. Equation (4.13) becomes a mere definition of the currents $I^Y$ ensuring that $[Y]$ follow the chemostat-imposed protocol.

*Example.* The reaction rates of the CRN (4.3) (represented in Fig. 4.1) are given by

$$r_1^+([Z]) = k_1^+[E][F], \qquad (4.14)$$

$$r_1^-([Z]) = k_1^-[EF], \qquad (4.15)$$

$$r_2^+([Z]) = k_2^+[EF], \qquad (4.16)$$

$$r_2^-([Z]) = k_2^-[EW], \qquad (4.17)$$

$$r_3^+([Z]) = k_3^+[EW], \qquad (4.18)$$

$$r_3^-([Z]) = k_3^-[E][W], \qquad (4.19)$$

$$r_4^+([Z]) = k_4^+[EF][S], \qquad (4.20)$$

$$r_4^-([Z]) = k_4^-[EFS], \qquad (4.21)$$

$$r_5^+([Z]) = k_5^+[EFS], \qquad (4.22)$$

$$r_5^-([Z]) = k_5^-[EW][P]. \qquad (4.23)$$
<figure><img src="image_41.png" /><figcaption>Figure 4.1: Pictorial illustration of the CRN (4.3) when the species F, S, W, and P are exchanged with chemostats. Black (numbered) arrows represent chemical reactions, while green arrows crossing the gray boundary represent the exchange processes with the chemostats (represented by flasks). Here, all reactions are assumed to be reversible even if they are represented with a hypergraph notation by single arrows.</figcaption></figure>

If we now assume that the species P, W, S, and F are exchanged (like in Fig. 4.1), the exchange current vector reads

$$ \mathbf{I} = \begin{pmatrix} E & EF & EW & EFS & P & W & S & F \\ 0 & 0 & 0 & 0 & I_P & I_W & I_S & I_F \end{pmatrix} \quad (4.24) $$

### 4.3.2 Steady State

The steady state $[Z]_{ss}$ of the rate equation (4.7), defined as

$$ S\mathbf{J}([Z]_{ss}) + \mathbf{I} = 0, \quad (4.25) $$

is said to be an equilibrium steady state if

$$ \mathbf{J}([Z]_{\text{eq}}) = 0, \quad (4.26) $$

which implies that $\mathbf{I} = 0$.

*Remark.* When CRNs have a well defined equilibrium state, they are said to be *detailed balanced* (not to be confused with *local detailed balance* which will be discussed in Sec. 4.4).

### 4.3.3 Conservation Laws

The $n_\ell \le n_s$ (linearly independent) left-null eigenvectors of the stoichiometric matrix, identified by the index $\lambda$,

$$ \ell^{\lambda} \cdot S = 0, \quad (4.27) $$
are named conservation laws. Indeed, each scalar

$$L^{\lambda} = \ell^{\lambda} \cdot [Z] \qquad (4.28)$$

is a conserved quantity of the rate equation (4.7) if CRNs were closed, i.e., $I = 0$:

$$\frac{d}{dt} L^{\lambda} = \ell^{\lambda} \cdot \frac{d}{dt} [Z] = \ell^{\lambda} \cdot S J([Z]) = 0. \qquad (4.29)$$

*Physical Remark.* The conservation laws identify fragments of (or entire) molecules, named, *moieties*, that remain intact in all reactions. The corresponding conserved quantities quantify the concentration of these moieties.

*Mathematical Remark.* The set of conservation laws {$\ell^{\lambda}$} is not unique: a linear combination of conservation laws is still a left-null eigenvector of the stoichiometric matrix. Different sets identify different moieties whose physical interpretation might not be obvious.

*Example.* The stoichiometric matrix (4.4) of the CRN (4.3) admits the following three conservation laws

$$\ell^{\text{E}} = \begin{pmatrix} E & 1 \\ EF & 1 \\ EW & 1 \\ EFS & 1 \\ P & 0 \\ W & 0 \\ S & 0 \\ F & 0 \end{pmatrix}, \quad \ell^{\text{S}} = \begin{pmatrix} E & 0 \\ EF & 0 \\ EW & 0 \\ EFS & 1 \\ P & 1 \\ W & 0 \\ S & 1 \\ F & 0 \end{pmatrix}, \quad \ell^{\text{F}} = \begin{pmatrix} E & 0 \\ EF & 1 \\ EW & 1 \\ EFS & 1 \\ P & 0 \\ W & 1 \\ S & 0 \\ F & 1 \end{pmatrix}. \qquad (4.30)$$

The concentrations of the corresponding moieties, i.e., the enzyme, the substrate and the fuel moieties, are given by $L^{\text{E}} = [\text{E}] + [\text{EF}] + [\text{EW}] + [\text{EFS}]$, $L^{\text{S}} = [\text{EFS}] + [\text{P}] + [\text{S}]$, $L^{\text{F}} = [\text{EF}] + [\text{EW}] + [\text{EFS}] + [\text{W}] + [\text{F}]$. The fragments corresponding to these moieties are highlighted in Fig. 4.2 with different colors.

## 4.4 Local Detailed Balance

The local detailed balance condition establishes a correspondence between the dynamics and the thermodynamics, and is crucial to build a nonequilibrium thermodynamic theory on top of a dynamical system. We introduce it here by comparing the dynamic and thermodynamic equilibrium conditions. We then generalize it to nonequilibrium conditions.

*Remark.* The general idea of our derivation of the local detailed balance is the following: i) Dynamic equilibrium relates $[Z]_{\text{eq}}$ to {$k_{\rho}^{\pm}$}; ii) Thermodynamic equilibrium relates $[Z]_{\text{eq}}$ to the standard chemical potentials {$\mu_{\alpha}^{\circ}(T)$}; iii) Local detailed balance relates {$k_{\rho}^{\pm}$} (dynamics) to {$\mu_{\alpha}^{\circ}(T)$} (thermodynamics); iv) Generalization to nonequilibrium.

### 4.4.1 Dynamic equilibrium

The steady state of closed CRNs must be an equilibrium steady state for thermodynamic consistency: if we do not provide (free) energy to dissipative systems such as CRNs, they must reach an equilibrium steady state.

By combining Eqs. (4.8), (4.9) and (4.26), we obtain

$$\frac{k_{\rho}^{+}}{k_{\rho}^{-}} = [Z]_{\text{eq}}^{S_{\rho}}, \qquad (4.31)$$
<figure><img src="image_42.png" /><figcaption>Figure 4.2: Pictorial illustration of the CRN (4.3) where the different moieties identified by the conservation laws in Eq. (4.30) are highlighted with different colors: orange for $\ell^E$, purple for $\ell^S$, and blue for $\ell^F$.</figcaption></figure>

for every $\rho \in \mathcal{R}$ or, equivalently,

$$RT \log \frac{k_{\rho}^{+}}{k_{\rho}^{-}} = RT \log([\mathbf{Z}]_{\text{eq}}) \cdot \mathbf{S}_{\rho}. \qquad (4.32)$$

Note that hereafter $\log \mathbf{v} = (\dots, \log v_i, \dots)^{\mathrm{T}}$ for very vector $\mathbf{v} = (\dots, v_i, \dots)^{\mathrm{T}}$.

### 4.4.2 Thermodynamic equilibrium

We consider here systems with temperature and pressure fixed by the environment they are exposed to. Equilibrium thermodynamics then dictates that the equilibrium state is a global minimum of the Gibbs free energy (density) that, for ideal mixture of chemical species, reads [72]

$$G([\mathbf{Z}]_{\text{eq}}, [Z_0]) = \sum_{\alpha \in \mathcal{S}} \mu_{\alpha}([\mathbf{Z}]_{\text{eq}}, [Z_0]) [Z_{\alpha}]_{\text{eq}} + \mu_0([\mathbf{Z}]_{\text{eq}}, [Z_0]) [Z_0] \quad (4.33)$$

with

$$\mu_{\alpha}([\mathbf{Z}]_{\text{eq}}, [Z_0]) = \hat{\mu}_{\alpha}^{\circ}(T) + RT \log \frac{[Z_{\alpha}]_{\text{eq}}}{\sum_{\beta \in \mathcal{S}} [Z_{\beta}]_{\text{eq}} + [Z_0]}, \quad (4.34)$$

$$\mu_0([\mathbf{Z}]_{\text{eq}}, [Z_0]) = \hat{\mu}_0^{\circ}(T) + RT \log \frac{[Z_0]}{\sum_{\beta \in \mathcal{S}} [Z_{\beta}]_{\text{eq}} + [Z_0]}, \quad (4.35)$$

being the chemical potential of $Z_{\alpha}$ and of the solvent, respectively; $\hat{\mu}_{\alpha}^{\circ}(T)$ and $\hat{\mu}_{0}^{\circ}(T)$ being the corresponding temperature-dependent standard chemical potentials; and $R$ and $T$ being the gas constant and the temperature, respectively.
For consistency with mass-action kinetics, we consider the dilute solution limit, i.e., $[Z_0] \gg \sum_{\alpha \in S}[Z_\alpha]$, which allows us to approximate

$$ \log \left( \frac{[Z_\alpha]_{\text{eq}}}{\sum_{\beta \in S} [Z_\beta]_{\text{eq}} + [Z_0]} \right) [Z_\alpha]_{\text{eq}} \approx \log \left( \frac{[Z_\alpha]_{\text{eq}}}{[Z_0]} \right) [Z_\alpha]_{\text{eq}}, \quad (4.36) $$

$$ \log \left( \frac{[Z_0]}{\sum_{\beta \in S} [Z_\beta]_{\text{eq}} + [Z_0]} \right) [Z_0] \approx - \sum_{\beta \in S} [Z_\beta]_{\text{eq}}, \quad (4.37) $$

using a Taylor expansion $^{10}$. This implies that the chemical potential of $Z_\alpha$ depends only on its concentration

$$ \mu_{\alpha}([Z_{\alpha}]_{\text{eq}}) \approx \mu_{\alpha}^{\circ}(T) + RT \log[Z_{\alpha}]_{\text{eq}}, \quad (4.39) $$

by rescaling the standard chemical potentials according to

$$ \mu_{\alpha}^{\circ}(T) := \hat{\mu}_{\alpha}^{\circ}(T) - RT \log[Z_{0}], \quad (4.40) $$

since $[Z_0]$ is a constant quantity. Hence, for ideal dilute solutions, the free energy reads

$$ G([Z]_{\text{eq}}) = \sum_{\alpha \in S} (\mu_{\alpha}([Z_{\alpha}]_{\text{eq}}) - RT)[Z_{\alpha}]_{\text{eq}} + \hat{\mu}_{0}^{\circ}(T)[Z_{0}]. \quad (4.41) $$

*Mathematical Remark.* Note that i) the last term in Eq. (4.41), i.e., $\hat{\mu}_{0}^{\circ}(T)[Z_0]$, represents a constant shift of the Gibbs free energy and will be neglected hereafter, and ii)

$$ \mu_{\alpha}([Z_{\alpha}]_{\text{eq}}) = \frac{\partial G([Z]_{\text{eq}})}{\partial [Z_{\alpha}]_{\text{eq}}}. \quad (4.42) $$

Equilibrium thermodynamics dictates that $[Z]_{\text{eq}}$ is a global minimum of $G$. Thus,

$$ \lim_{\epsilon \to 0} \frac{G([Z]_{\text{eq}} + \epsilon S_{\rho}) - G([Z]_{\text{eq}})}{\epsilon} = 0, \quad (4.43) $$

for every $\rho$ or equivalently,

$$ \Delta_{\rho} G([Z]_{\text{eq}}) := \mu([Z]_{\text{eq}}) \cdot S_{\rho} = 0, \quad (4.44) $$

where $\mu([Z]_{\text{eq}}) = (\dots, \mu_{\alpha}([Z_{\alpha}]_{\text{eq}}), \dots)_{\alpha \in S}^{\text{T}}$. Equation (4.44) can be rewritten as

$$ -\mu^{\circ}(T) \cdot S_{\rho} = RT \log([Z]_{\text{eq}}) \cdot S_{\rho}, \quad (4.45) $$

with $\mu^{\circ}(T) = (\dots, \mu_{\alpha}^{\circ}(T), \dots)_{\alpha \in S}^{\text{T}}$.

*Mathematical Remark.* Equation (4.43) accounts for all possible changes of the concentrations that are consistent with the conservation laws. An alternative derivation of Eq. (4.44) uses Lagrange multipliers $^{11}$.

$$ \log \left( \frac{[Z_0]}{\epsilon + [Z_0]} \right) = \epsilon \left[ \frac{\epsilon + [Z_0]}{[Z_0]} \frac{(-[Z_0])}{(\epsilon + [Z_0])^2} \right]_{\epsilon=0} + O(\epsilon^2), \quad (4.38) $$

truncated at the first order.

$$ \mathcal{L}([Z]_{\text{eq}}, f) = G([Z]_{\text{eq}}) - \sum_{\lambda} f_{\lambda} C^{\lambda}([Z]_{\text{eq}}), \quad (4.46) $$

<footnote>$^{10}$By defining $\epsilon = \sum_{\beta \in S} [Z_{\beta}]_{\text{eq}}$, we use the Taylor expansion</footnote>

<footnote>$^{11}$Given the Lagrangian function</footnote>
### 4.4.3 Equilibrium Local Detailed Balance

By combining Eqs. (4.32) and (4.45), we obtain the local detailed balance condition

$$RT \log \frac{k_{\rho}^{+}}{k_{\rho}^{-}} = -\boldsymbol{\mu}^{\circ}(T) \cdot \boldsymbol{S}_{\rho}. \qquad (4.51)$$

*Physical Remark.* Our derivation of the local detailed balance is only based on the assumption of thermodynamic consistency, namely, the rate equation (4.7) for closed CRNs must admit an equilibrium steady state $[Z]_{eq}$ which corresponds to the global minimum of the (equilibrium) Gibbs free energy (4.41).

*Remark.* Note the difference between detailed balanced CRNs discussed in Subs. 4.3.2, namely, the existence of a well defined equilibrium steady state, and the local detailed balance condition derived here, namely, the correspondence between the kinetic constants and the standard chemical potentials.

### 4.4.4 Non-Equilibrium Local Detailed Balance

We now generalize the local detailed balance to nonequilibrium conditions. To do so, we assume that all degrees of freedom other than concentrations are always equilibrated. The temperature $T$ is set by the solvent playing the role of a thermal reservoir, the pressure $p$ is set by the environment the system is exposed to, and diffusion processes are fast enough to keep the chemical species homogeneously distributed. This physically means that CRNs would be an equilibrated mixture in absence of reactions for every concentration $[Z]$. The dynamics is merely interconverting the mixture $[Z](t)$ into another $[Z](t+dt)$. In this way, thermodynamic state functions, namely, the Gibbs free energy and the chemical potentials, can be specified by their equilibrium form (in Eqs. (4.41) and (4.39), respectively) but expressed in terms of nonequilibrium concentrations $[Z]$:

$$G([Z]) = \sum_{\alpha \in S} (\mu_{\alpha}([Z_{\alpha}]) - RT)[Z_{\alpha}], \qquad (4.52)$$

with

$$\mu_{\alpha}([Z_{\alpha}]) = \mu_{\alpha}^{\circ}(T) + RT \log[Z_{\alpha}]. \qquad (4.53)$$

We analogously introduce

$$\Delta_{\rho}G([Z]) := \mu([Z]) \cdot S_{\rho} \neq 0, \qquad (4.54)$$

which is in general different from zero in nonequilibrium conditions.

where

$$C^{\lambda}([Z]_{\text{eq}}) = \ell^{\lambda} \cdot [Z]_{\text{eq}} - \bar{L}^{\lambda}, \qquad (4.47)$$

with $\{\bar{L}^{\lambda}\}$ being the actual values of the conserved quantities and $f = (..., f_{\lambda},...)$ being the Lagrange multipliers, the global minimum of $G$ consistent with the conservation laws satisfies

$$\frac{\partial \mathcal{L}([Z]_{\text{eq}}, f)}{\partial f_{\lambda}} = -C^{\lambda}([Z]_{\text{eq}}) = 0, \qquad (4.48)$$

$$\frac{\partial \mathcal{L}([Z]_{\text{eq}}, f)}{\partial [Z_{\alpha}]_{\text{eq}}} = \mu_{\alpha}([Z_{\alpha}]_{\text{eq}}) - \sum_{\lambda} f_{\lambda} \ell_{\alpha}^{\lambda} = 0, \qquad (4.49)$$

with $\ell_{\alpha}^{\lambda}$ being the $\alpha$ entry of $\ell^{\lambda}$. The condition in Eq. (4.49) can be rewritten as

$$\mu([Z]_{\text{eq}}) = \sum_{\lambda} f_{\lambda} \ell^{\lambda}, \qquad (4.50)$$

which implies Eq. (4.44) because of Eq. (4.27).
Thus, by summing $-RT \log[\mathbf{Z}] \cdot \mathbf{S}_{\rho}$ (with $[\mathbf{Z}]$ the nonequilibrium concentrations) on both sides of Eq. (4.51), and using Eqs. (4.9) and (4.54), we can formulate the local detailed balance in terms of a flux-force relation:

$$RT \log \frac{r_{\rho}^{+}([\mathbf{Z}])}{r_{\rho}^{-}([\mathbf{Z}])} = -\Delta_{\rho}G([\mathbf{Z}]), \qquad (4.55)$$

where $\Delta_{\rho}G$ plays the role of the thermodynamic force driving the reaction. Indeed, if $\Delta_{\rho}G < 0$, then $r_{\rho}^{+} > r_{\rho}^{-}$ and reaction $\rho$ proceeds in the forward direction $J_{\rho} > 0$. On the other hand, if $\Delta_{\rho}G > 0$, then $r_{\rho}^{+} < r_{\rho}^{-}$ and reaction $\rho$ proceeds in the backward direction $J_{\rho} < 0$.

## 4.5 Thermodynamics of Closed CRNs

We now formulate the first and the second law of thermodynamics for closed CRNs. We will use the second law and the conservation laws to prove that the concentrations will eventually reach equilibrium, i.e., $\lim_{t \to +\infty} [\mathbf{Z}](t) = [\mathbf{Z}]_{\text{eq}}$, for every initial condition $[\mathbf{Z}](0)$.

### 4.5.1 First Law

We introduce the enthalpy $H$, quantifying the “internal” energy of CRNs$^{12}$, as

$$H([\mathbf{Z}]) := \frac{\partial(G([\mathbf{Z}])/T)}{\partial 1/T} = \mathbf{h}^{\circ}(T) \cdot [\mathbf{Z}], \qquad (4.57)$$

where $\mathbf{h}^{\circ}(T) = (\dots, h_{\alpha}^{\circ}(T), \dots)_{\alpha \in S}^{\mathrm{T}}$ and

$$h_{\alpha}^{\circ}(T) := \frac{\partial(\mu_{\alpha}^{\circ}(T)/T)}{(\partial 1/T)} = \frac{\partial(\mu_{\alpha}([Z_{\alpha}])/T)}{(\partial 1/T)} \qquad (4.58)$$

is the standard molar enthalpy of $Z_{\alpha}$. By taking the time derivative of $H([\mathbf{Z}])$ according to the rate equation (4.7), we obtain the nonequilibrium formulation of the first law of thermodynamics for closed CRNs:

$$\frac{d}{dt} H([\mathbf{Z}]) = \mathbf{h}^{\circ}(T) \cdot \mathcal{S} \mathbf{J}([\mathbf{Z}]) =: \dot{Q}. \qquad (4.59)$$

Since the heat exchange with the thermal reservoir is the only mechanism to exchange energy, we recognize the heat flux $\dot{Q}$ on the rightmost-hand side of Eq. (4.59).

### 4.5.2 Second Law

We introduce the entropy $S$ as

$$S([\mathbf{Z}]) := -\frac{\partial G([\mathbf{Z}])}{\partial T} = \sum_{\alpha \in S} (s_{\alpha}([Z_{\alpha}]) + R)[Z_{\alpha}], \qquad (4.60)$$

with

$$s_{\alpha}([Z_{\alpha}]) := s_{\alpha}^{\circ}(T) - R \log[Z_{\alpha}], \qquad (4.61)$$

$$s_{\alpha}^{\circ}(T) := -\partial \mu_{\alpha}^{\circ}(T)/\partial T, \qquad (4.62)$$

$$H(S, p, [\mathbf{Z}]) := [U(S, V, [\mathbf{Z}]) + pV]_{V=V(S, p, [\mathbf{Z}])}, \qquad (4.56)$$

with $V(S, p, [\mathbf{Z}])$ the function expressing the volume in terms of entropy $S$, pressure $p$, and concentrations $[\mathbf{Z}]$.

<footnote>$^{12}$The enthalpy $H$ is actually the Legendre transform of the internal energy $U$. The natural variables of $U$ are entropy $S$, volume $V$, and concentrations $[\mathbf{Z}]$, and enthalpy is defined as</footnote>
being the molar entropy and the standard molar entropy of $Z_\alpha$, respectively. The time derivative of $S$ according to the rate equation (4.7) reads

$$ \frac{d}{dt} S([\mathbf{Z}]) = \sum_{\alpha \in S} s_{\alpha}([Z_{\alpha}]) \frac{d}{dt} [Z_{\alpha}] \quad (4.63a) $$

$$ = s([\mathbf{Z}]) \cdot \mathcal{S}\mathbf{J}([\mathbf{Z}]) \quad (4.63b) $$

$$ = T^{-1}(h^{\circ}(T) - \mu([\mathbf{Z}])) \cdot \mathcal{S}\mathbf{J}([\mathbf{Z}]) \quad (4.63c) $$

with $s([\mathbf{Z}]) = (\dots, s_\alpha([Z_\alpha]), \dots)^T_{\alpha \in S}$ and $s([\mathbf{Z}]) = T^{-1}(h^{\circ}(T) - \mu([\mathbf{Z}]))$ $^{13}$. By recognizing the (reversible) entropy flow

$$ \dot{S}_e := T^{-1}h^{\circ}(T) \cdot \mathcal{S}\mathbf{J}([\mathbf{Z}]) = \frac{\dot{Q}}{T}, \quad (4.65) $$

representing the (reversible) entropy changes in the environment, and the entropy production rate

$$ \dot{\Sigma} := -T^{-1}\boldsymbol{\mu}([\mathbf{Z}]) \cdot \mathcal{S}\mathbf{J}([\mathbf{Z}]) \geq 0, \quad (4.66) $$

we obtain the nonequilibrium formulation of the second law of thermodynamics for closed CRNs:

$$ \frac{d}{dt} S([\mathbf{Z}]) = \dot{S}_e + \dot{\Sigma}. \quad (4.67) $$

Note that, by using Eq. (4.54), the entropy production rate can be written as

$$ \dot{\Sigma} = -T^{-1} \sum_{\rho \in \mathcal{R}} \Delta_{\rho} G([\mathbf{Z}]) J_{\rho}([\mathbf{Z}]) \geq 0, \quad (4.68) $$

which, by using the local detailed balance condition (4.55), can be rewritten in a manifestly non-negative form

$$ \dot{\Sigma} = R \sum_{\rho \in \mathcal{R}} \underbrace{(r_{\rho}^{+}([\mathbf{Z}]) - r_{\rho}^{-}([\mathbf{Z}])) \log \frac{r_{\rho}^{+}([\mathbf{Z}])}{r_{\rho}^{-}([\mathbf{Z}])}}_{=-\Delta_{\rho}G([\mathbf{Z}])J_{\rho}([\mathbf{Z}])/(RT)} \geq 0. \quad (4.69) $$

We so notice that each term $-\Delta_{\rho}G([\mathbf{Z}])J_{\rho}([\mathbf{Z}])$ of the sum in Eq. (4.68) is greater than or equal to zero: the thermodynamic forces and currents are always aligned, namely, $-\Delta_{\rho}G > 0$ (resp. $-\Delta_{\rho}G < 0$) if and only if $J_{\rho} > 0$ (resp. $J_{\rho} < 0$).

*Remark.* The entropy production rate vanishes only at equilibrium, namely, when $\Delta_{\rho}G([\mathbf{Z}]_{\text{eq}}) = 0$ and $J_{\rho}([\mathbf{Z}]_{\text{eq}}) = 0$ for every reaction $\rho$.

$$ h_{\alpha}^{\circ}(T) = \frac{1}{T} \frac{\partial \mu_{\alpha}}{\partial 1/T} + \mu_{\alpha} \frac{\partial 1/T}{\partial 1/T} \quad (4.64a) $$

$$ = \frac{1}{T} \frac{\partial \mu_{\alpha}}{\partial T} \frac{\partial T}{\partial 1/T} + \mu_{\alpha} \frac{\partial 1/T}{\partial 1/T} \quad (4.64b) $$

$$ = \frac{1}{T} s_{\alpha}(-T^2) + \mu_{\alpha} \quad (4.64c) $$

proving that $s([\mathbf{Z}]) = T^{-1}(h^{\circ}(T) - \mu([\mathbf{Z}]))$.

<footnote>$^{13}$Indeed, $h^{\circ}(T) = \partial(\mu_{\alpha}/T)/(\partial 1/T)$ implies</footnote>
### 4.5.3 Thermodynamic Potential

We now show that the Gibbs free energy in Eq. (4.52) is the proper thermodynamic potential of closed CRNs, namely, it monotonically decreases during the dynamics and it is lower bounded by its equilibrium value given in Eq. (4.41).

First, we examine the time derivative of $G([Z])$ and we obtain

$$ \frac{d}{dt}G([Z]) = \sum_{\alpha \in S} \mu_{\alpha}([Z_{\alpha}]) \frac{d}{dt}[Z_{\alpha}] \quad (4.70a) $$

$$ = \mu([Z]) \cdot \mathcal{S}\mathcal{J}([Z]) \quad (4.70b) $$

and, by using the definition of the entropy production rate in Eq. (4.66),

$$ \frac{d}{dt}G([Z]) = -T\dot{\Sigma} \le 0, \quad (4.71) $$

which physically means that the dynamics dissipates the Gibbs free energy.

*Remark.* Equation (4.71) is an alternative formulation of the second law (4.67) for closed CRNs.

Second, we examine the difference $G([Z]) - G([Z]_{eq})$ (with $G([Z])$ and $G([Z]_{eq})$ given in Eqs. (4.52) and (4.41), respectively). To do so, we start by inspecting Eq. (4.44). We notice that the vector of the equilibrium chemical potentials belongs to the cokernel of the stoichiometric matrix. Hence, it can be written as a linear combination of conservation laws (defined in Eq. (4.27)):

$$ \mu([Z]_{eq}) = \sum_{\lambda} f_{\lambda} \ell^{\lambda}, \quad (4.72) $$

with {$f_{\lambda}$} some coefficients. Consequently,

$$ \mu([Z]_{eq}) \cdot [Z]_{eq} = \sum_{\lambda} f_{\lambda} \underbrace{\ell^{\lambda} \cdot [Z]_{eq}}_{=L^{\lambda}} \quad (4.73a) $$

$$ = \sum_{\lambda} f_{\lambda} \underbrace{\ell^{\lambda} \cdot [Z]}_{=L^{\lambda}} \quad (4.73b) $$

$$ = \mu([Z]_{eq}) \cdot [Z]. \quad (4.73c) $$

This implies that

$$ \mu([Z]) \cdot [Z] - \mu([Z]_{eq}) \cdot [Z]_{eq} = RT \sum_{\alpha \in S} [Z_{\alpha}] \log \frac{[Z_{\alpha}]}{[Z_{\alpha}]_{eq}}, \quad (4.74) $$

and we can write the difference $G([Z]) - G([Z]_{eq})$ as a relative entropy (or Kullback-Leibler divergence) for non-normalized distributions:

$$ G([Z]) - G([Z]_{eq}) = RT\mathcal{D}([Z] \| [Z]_{eq}) \ge 0, \quad (4.75) $$

with

$$ \mathcal{D}([Z] \| [Z]_{eq}) := \sum_{\alpha \in S} \left( [Z_{\alpha}] \log \frac{[Z_{\alpha}]}{[Z_{\alpha}]_{eq}} - ([Z_{\alpha}] - [Z_{\alpha}]_{eq}) \right) \ge 0. \quad (4.76) $$

Note that $\mathcal{D}([Z] \| [Z]_{eq}) \ge 0$ because of the logarithmic inequality $\log x \le x - 1$ (when $x = [Z_{\alpha}]/[Z_{\alpha}]_{eq}$).

*Mathematical Remark.* Equations (4.71) and (4.75) prove that the Gibbs free energy plays the role of a Lyapunov function for the closed system dynamics.
## 4.6 Thermodynamics of Open CRNs

We now generalize the results obtained in Sec. 4.5 to open CRNs. Before doing that, we first introduce the thermodynamic meaning of the exchange currents entering the rate equation (4.7).

### 4.6.1 Chemostats

In these lecture notes, we consider CRNs coupled to chemostats. This might be interpreted as if each exchanged species $Z_\alpha$ is involved in a chemical reaction (labeled $\rho_\alpha^{ex}$) like

$$Y_{\alpha} \underset{\alpha \in S_Y}{\overset{\rho_{\alpha}^{\text{ex}}}{\rightleftharpoons}} Y_{\alpha}^{\text{ex}}, \qquad (4.77)$$

with the corresponding chemostats $Y_\alpha^{ex}$. These reactions are assumed to be always equilibrated. This thermodynamically means that the chemical potential of each exchanged species is controlled by the corresponding chemostats:

$$\mu_{\alpha}([Y_{\alpha}]) = \mu_{\alpha}^{\text{ex}}, \qquad (4.78)$$

for every $\alpha \in S_Y$.

*Remark.* The chemical potentials of the exchanged species $\mu_Y = (..., \mu_\alpha, ...)^T_{\alpha \in S_Y}$ are thus externally controlled parameters like their concentrations $[Y]$.

### 4.6.2 First and Second Law

When the chemical species are exchanged with the environment, by using the rate equation (4.7), the first law (4.59) becomes

$$\frac{d}{dt} H([Z]) = \dot{Q} + \sum_{\alpha \in S_Y} h_{\alpha}^{\circ}(T) I_{\alpha}. \qquad (4.79)$$

Therefore, while the expression of the second law (4.67) remains the same, the entropy flow (4.65) becomes

$$\dot{S}_e = \frac{\dot{Q}}{T} + \sum_{\alpha \in S_Y} s_{\alpha}([Z_{\alpha}]) I_{\alpha}. \qquad (4.80)$$

The second term on the right-hand side of Eq. (4.79) (resp. Eq. (4.80)) accounts for the molar energy (resp. entropy) exchanged with the environment through the exchange of chemical species.

The formulation of the second law in Eq. (4.71) becomes

$$\frac{d}{dt} G([Z]) = -T \dot{\Sigma} + \sum_{\alpha \in S_Y} \mu_{\alpha}([Z_{\alpha}]) I_{\alpha}, \qquad (4.81)$$

where

$$\dot{w}_{\text{chem}} := \sum_{\alpha \in S_Y} \mu_{\alpha}([Z_{\alpha}]) I_{\alpha}. \qquad (4.82)$$

is the chemical work rate quantifying the free energy exchanged with the environment through the exchange of chemical species. By providing free energy, the chemical work can balance dissipation and maintain CRNs out of equilibrium.

*Physical Remark.* These formulations of the first and second laws do not provide any precise guideline to understand which mechanisms (if any) maintain CRNs out of equilibrium (besides the exchange processes with the chemostats). To identify these mechanisms, we need to use the conservation laws.
### 4.6.3 Broken Conservation Laws

In open CRNs, some of the quantities defined in Eq. (4.28) are not conserved anymore. The corresponding conservation laws are said to be broken. The other conservation laws, which still correspond to conserved quantities, are said to be unbroken. To systematically identify broken and unbroken conservation laws, we proceed as follows.

First, we consider the $n_u \le n_X$ (lineally independent) left-null eigenvectors of the (sub)stoichiometric matrix $S^X$, identified by the index $\lambda_u$:

$$ \ell_X^{\lambda_u} \cdot S^X = 0, \qquad (4.83) $$

where the subscript $X$ stresses that $\ell_X^{\lambda_u}$ has $n_X$ entries (corresponding to the internal species). The vectors

$$ \ell^{\lambda_u} := \begin{pmatrix} \ell_X^{\lambda_u} \\ \mathbf{0}_Y \end{pmatrix}, \qquad (4.84) $$

with $\mathbf{0}_Y$ the vector with $n_Y$ null entries, are unbroken conservation laws. Indeed,

$$ \ell^{\lambda_u} \cdot S = \ell_X^{\lambda_u} \cdot S^X + \mathbf{0}_Y \cdot S^Y = 0, \qquad (4.85) $$

because of Eq. (4.83), and the quantities

$$ L^{\lambda_u} = \ell^{\lambda_u} \cdot [\mathbf{Z}] \qquad (4.86) $$

are conserved quantities in open CRNs:

$$ \frac{d}{dt} L^{\lambda_u} = \ell^{\lambda_u} \cdot \frac{d}{dt} [\mathbf{Z}] = \underbrace{\ell^{\lambda_u} \cdot S \mathbf{J}([\mathbf{Z}])}_{=0} + \underbrace{\mathbf{0}_Y \cdot \mathbf{I}^Y}_{=0} = 0, \qquad (4.87) $$

where we used Eq. (4.7).

Second, we use the set $\{\ell^{\lambda_u}\}$ of unbroken conservation laws to express the set of conservation laws $\{\ell^{\lambda}\}$ according to

$$ \{\ell^{\lambda}\} = \{\ell^{\lambda_u}\} \cup \{\ell^{\lambda_b}\} \qquad (4.88) $$

where $\{\ell^{\lambda_b}\}$ is the set of $n_b$ broken conservation laws (with $n_\ell = n_u + n_b$). Indeed, each vector $\ell_Y^{\lambda_b} = (\dots, \ell_\alpha^{\lambda_b}, \dots)^T_{\alpha \in S_Y}$ (collecting the entries of $\ell^{\lambda_b}$ corresponding to the exchanged species) satisfies $\ell_Y^{\lambda_b} \neq \mathbf{0}_Y$, otherwise $\ell^{\lambda_b}$ would be an unbroken conservation law. Hence, the quantities

$$ L^{\lambda_b} = \ell^{\lambda_b} \cdot [\mathbf{Z}] \qquad (4.89) $$

are in general not conserved (and their variation is solely due to the exchanges with the chemostats):

$$ \frac{d}{dt} L^{\lambda_b} = \ell^{\lambda_b} \cdot \frac{d}{dt} [\mathbf{Z}] = \underbrace{\ell^{\lambda_b} \cdot S \mathbf{J}([\mathbf{Z}])}_{=0} + \underbrace{\ell_Y^{\lambda_b} \cdot \mathbf{I}^Y}_{\neq 0} \neq 0. \qquad (4.90) $$

*Example.* When the species P, W, S, and F of the CRN (4.3) are exchanged, the conservation laws $\ell^S$ and $\ell^F$ in Eq. (4.30) are broken, while the conservation law $\ell^E$ in Eq. (4.30) is unbroken. In Fig. 4.3, we show that the moieties corresponding to the conservation laws $\ell^S$ and $\ell^F$ are exchanged with chemostats.
<figure><img src="image_43.png" /><figcaption>Figure 4.3: Pictorial illustration of the CRN (4.3) where the moieties corresponding to the conservation laws $\ell^S$ and $\ell^F$ in Eq. (4.30) are exchanged with chemostats. Crucially, the moiety $\ell^S$ (resp. $\ell^F$) is exchanged with the purple (resp. blue) chemostats exchanging S and P (resp. F and W).</figcaption></figure>

### 4.6.4 Force and Potential Exchanged Species

In open CRNs the same moiety (the fragment of molecules identified by a conservation law) can be exchanged with more than one chemostat (as shown in Fig. 4.3), which can be identified by using the broken conservation laws.

We split the set of exchanged species $\{Y_\alpha\}$, and the corresponding set of indexes $S_Y$, into two disjoint subsets: the set of $n_{Y_p}$ potential species $\{Y_\alpha\}$ with $\alpha \in S_{Y_p}$, and the set of $n_{Y_f}$ force species $\{Y_\alpha\}$ with $\alpha \in S_{Y_f}$. The former is defined as the smallest subset of exchanged species such that all the vectors $\{\ell^{\lambda_b}\}$ — independently of the specific representation — satisfy

$$ \ell_{Y_p}^{\lambda_b} \neq 0_{Y_p}, \qquad (4.91) $$

with $\ell_{Y_p}^{\lambda_b} = (\dots, \ell_{\alpha}^{\lambda_b}, \dots)^T_{\alpha \in S_{Y_p}}$ (collecting the entries of $\ell^{\lambda_b}$ corresponding to the potential species), and $0_{Y_p}$ being the vector with $n_{Y_p}$ null entries.

**Mathematical Remark.** The vectors $\{\ell_{Y_p}^{\lambda_b}\}$ are linearly independent. Indeed, if they were linearly dependent, there would a representation of $\{\ell^{\lambda_b}\}$ such that

$$ \ell_{Y_p}^{\lambda_b} = 0_{Y_p}, \qquad (4.92) $$

for at least one $\lambda_b$. This, together with $S_{Y_p}$ being the smallest subset such that Eq. (4.91) holds, implies that the number of broken conservations is equal to the number of potential species, i.e., $n_{Y_p} = n_b$.

**Mathematical Remark.** By applying the same splitting to the $(n_b \times n_s)$ matrix $\mathbb{L}^b$, collecting all the broken conservation laws as row vectors, we get

$$ \mathbb{L}^b = (\mathbb{L}_X^b, \mathbb{L}_{Y_f}^b, \mathbb{L}_{Y_p}^b), \qquad (4.93) $$

where $\mathbb{L}_{Y_p}^b$ is a $(n_b \times n_{Y_p})$ square and invertible matrix (since the vectors $\{\ell_{Y_p}^{\lambda_b}\}$ are linearly independent).
*Mathematical Remark.* By applying the same splitting to the stoichiometric matrix, we get

$$ S = \begin{pmatrix} S^X \\ S^{Y_f} \\ S^{Y_p} \end{pmatrix}. \qquad (4.94) $$

By removing all rows of the (sub)stoichiometric matrix $S^{Y_p}$ from $S$, the only linearly-independent left null eigenvectors of the resulting matrix satisfy <fnref n="14" />:

$$ \begin{pmatrix} \ell_X^{\lambda_u} \\ \mathbf{0}_{Y_f} \end{pmatrix}, \qquad (4.95) $$

with $\mathbf{0}_{Y_f}$ being the vector with $n_{Y_f}$ null entries.

*Physical Remark.* Every time a potential species is exchanged with a chemostat, a new conservation law is broken and a new moiety is exchanged with that chemostat. When a force species is exchanged with a chemostat, no new conservation laws are broken and that chemostat exchanges a moiety with the CRN that was already exchanged with another chemostat, thus establishing a moiety-flux between two (or more) chemostats.

*Remark.* The splitting of the exchanged species in force and potential species is not unique. Different choices have different physical interpretation.

*Example.* Possible sets of potential species of the CRN (4.3) when the species P, W, S, and F are exchanged are {P, W} implying

$$ \ell_{Y_p}^S = \begin{pmatrix} 1 \\ 0 \end{pmatrix}_W \quad \text{and} \quad \ell_{Y_p}^F = \begin{pmatrix} 0 \\ 1 \end{pmatrix}_W, \qquad (4.96) $$

or {P, F} implying

$$ \ell_{Y_p}^S = \begin{pmatrix} 1 \\ 0 \end{pmatrix}_F \quad \text{and} \quad \ell_{Y_p}^F = \begin{pmatrix} 0 \\ 1 \end{pmatrix}_F, \qquad (4.97) $$

or {S, W} implying

$$ \ell_{Y_p}^S = \begin{pmatrix} 0 \\ 1 \end{pmatrix}_S \quad \text{and} \quad \ell_{Y_p}^F = \begin{pmatrix} 1 \\ 0 \end{pmatrix}_S, \qquad (4.98) $$

or {S, F} implying

$$ \ell_{Y_p}^S = \begin{pmatrix} 1 \\ 0 \end{pmatrix}_F \quad \text{and} \quad \ell_{Y_p}^F = \begin{pmatrix} 0 \\ 1 \end{pmatrix}_F. \qquad (4.99) $$

Let us consider the species {S, F} as potential species. This physically means that when the species S (resp. F) is exchanged the conservation law $\ell^S$ (resp. $\ell^F$) is broken. If afterwards, also the species P (resp. W) is exchanged, no new conservation law is broken, but the moiety corresponding to the conservation law $\ell^S$ (resp. $\ell^F$) is now exchanged with two different chemostats (see Fig. 4.3). Notice that the matrix $\mathbb{L}^b$ now reads

$$ \mathbb{L}^b = \begin{pmatrix} \ell^S & 0 & 0 & 0 & 1 & 1 & 0 \\ \ell^F & 0 & 1 & 1 & 1 & 0 & 1 \end{pmatrix} \qquad (4.100) $$

where the vertical lines split $\mathbb{L}^b$ into $\mathbb{L}_X^b$, $\mathbb{L}_{Y_f}^b$, and $\mathbb{L}_{Y_p}^b$.

<footnote n="14">¹⁴This can be proven by repeating the same reasoning we used to identify the unbroken conservation laws in Subs. 4.6.3. Every time we remove a row of $S^{Y_p}$ from $S$, we break one conservation law because of Eq. (4.91) and the dimension of the cokernel of resulting matrix thus decreases.</footnote>
### 4.6.5 Decomposition of the Entropy Production Rate

We now decompose the entropy production into a potential change and two work contributions. One work contribution accounts for the flows of moieties maintained across the CRN breaking the detailed balance condition. The other is due to time-dependent changes in the externally controlled chemostated concentrations. This decomposition constitutes a new formulation of the second law for open CRNs.

We start by rewriting the entropy production rate in Eq. (4.66) by explicitly considering the contributions of the different species:

$$T\dot{\Sigma} := -(\boldsymbol{\mu}_X \cdot \mathbb{S}^X + \boldsymbol{\mu}_{Y_f} \cdot \mathbb{S}^{Y_f} + \boldsymbol{\mu}_{Y_p} \cdot \mathbb{S}^{Y_p}) \mathbf{J} \geq 0, \quad (4.101)$$

where $\boldsymbol{\mu}_i = (\dots, \mu_\alpha, \dots)_{\alpha \in i}$ with $i = \{S_X, S_{Y_f}, S_{Y_p}\}$, and we do not explicitly write the $[\mathbf{Z}]$ dependence of the chemical potentials and reaction currents for the sake of compactness.

*Remark.* Recall that the chemical potentials of the exchanged species, $\boldsymbol{\mu}_Y = (\boldsymbol{\mu}_{Y_f}, \boldsymbol{\mu}_{Y_p})^T$, are controlled by the chemostats.

We then need to recognize that part of the free energy externally exchanged through the currents $\mathbf{I}^Y$ modifies the free energy of the system, while another part is transferred through the system. To proceed, we exploit the conservation laws. We consider

$$\mathbb{L}^b \mathbb{S} = 0 = \mathbb{L}_X^b \mathbb{S}^X + \mathbb{L}_{Y_f}^b \mathbb{S}^{Y_f} + \mathbb{L}_{Y_p}^b \mathbb{S}^{Y_p}, \quad (4.102)$$

and using that $\mathbb{L}_{Y_p}^b$ can be inverted (as discussed in 4.6.4), we obtain

$$\mathbb{S}^{Y_p} = -(\mathbb{L}_{Y_p}^b)^{-1}(\mathbb{L}_X^b \mathbb{S}^X + \mathbb{L}_{Y_f}^b \mathbb{S}^{Y_f}), \quad (4.103)$$

binding the variation of the number of molecules of the potential species to the other species. By inserting Eq. (4.103) in Eq. (4.101), we get

$$T\dot{\Sigma} = -(\boldsymbol{\mu}_X \cdot \mathbb{S}^X + \boldsymbol{\mu}_{Y_f} \cdot \mathbb{S}^{Y_f} + \boldsymbol{\mu}_{Y_p} \cdot (\mathbb{L}_{Y_p}^b)^{-1}(\mathbb{L}_X^b \mathbb{S}^X + \mathbb{L}_{Y_f}^b \mathbb{S}^{Y_f}))\mathbf{J}. \quad (4.104)$$

By adding and subtracting $\boldsymbol{\mu}_{Y_p} \cdot (\mathbb{L}_{Y_p}^b)^{-1}\mathbb{L}_{Y_p}^b \mathbb{S}^{Y_p} = \boldsymbol{\mu}_{Y_p} \cdot \mathbb{S}^{Y_p}$, the entropy production rate becomes

$$T\dot{\Sigma} = -\mathcal{F} \cdot \mathbb{S}\mathbf{J}, \quad (4.105)$$

where

$$\mathcal{F} = (\boldsymbol{\mu} \cdot \mathbb{I} - \boldsymbol{\mu}_{Y_p} \cdot (\mathbb{L}_{Y_p}^b)^{-1}\mathbb{L}^b)^T \quad (4.106)$$

with $\mathbb{I}$ the identity matrix. Note that we just introduced a term in the expression of the entropy production rate, i.e., $\boldsymbol{\mu}_{Y_p} \cdot (\mathbb{L}_{Y_p}^b)^{-1}\mathbb{L}^b$, whose net contribution vanishes since $\mathbb{L}^b \mathbb{S} = 0$. Nevertheless, introducing this term allows us to recognize the proper thermodynamic potential of open CRNs.

To do so, we guess the following state function

$$\Psi([\mathbf{Z}]) = \mathcal{F} \cdot [\mathbf{Z}], \quad (4.107)$$

whose time derivative

$$\frac{d}{dt}\Psi([\mathbf{Z}]) = \mathcal{F} \cdot \frac{d}{dt}[\mathbf{Z}] + \frac{d}{dt}\mathcal{F} \cdot [\mathbf{Z}] \quad (4.108)$$

according to the rate equation (4.7) reads

$$\begin{aligned} \frac{d}{dt}\Psi([\mathbf{Z}]) = & \mathcal{F} \cdot \mathbb{S}\mathbf{J} + \mathcal{F}_Y \cdot \mathbf{I}^Y + RT \sum_{\alpha \in S} \frac{d}{dt}[Z_\alpha] \\ & - \frac{d}{dt}(\boldsymbol{\mu}_{Y_p}) \cdot (\mathbb{L}_{Y_p}^b)^{-1}\mathbb{L}^b[\mathbf{Z}] \end{aligned} \quad (4.109)$$
with $\mathcal{F}_Y = (\mu_Y \cdot \mathbb{1} - \mu_{Y_p} \cdot (\mathbb{1}_{Y_p}^b)^{-1}\mathbb{1}_{Y}^b)^T$ and $\mathbb{1}_{Y}^b = (\mathbb{1}_{Y_f}^b, \mathbb{1}_{Y_p}^b)$.

To proceed, i) we express $-\mathcal{F} \cdot \mathcal{S}\mathbf{J} = T\dot{\mathcal{S}}$ as a function of the other terms, ii) we recognize that $\mathcal{F}_Y \cdot \mathbf{I}^Y = \mathcal{F}_{Y_f} \cdot \mathbf{I}^{Y_f}$, and iii) we define $\mathcal{G}([\mathbf{Z}]) := \Psi([\mathbf{Z}]) - RT \sum_{\alpha \in \mathcal{S}} [Z_{\alpha}]$. We thus get

$$T\dot{\mathcal{S}} = -\frac{d}{dt}\mathcal{G}([\mathbf{Z}]) + \dot{w}_{\text{nc}} + \dot{w}_{\text{driv}} \geq 0. \qquad (4.110)$$

Here, we introduced the new free energy $\mathcal{G}$

$$\mathcal{G}([\mathbf{Z}]) = \mu \cdot [\mathbf{Z}] - RT \sum_{\alpha \in \mathcal{S}} [Z_{\alpha}] - \mu_{Y_p} \cdot (\mathbb{1}_{Y_p}^b)^{-1}\mathbb{1}^b[\mathbf{Z}] \qquad (4.111)$$

which can be written using the Gibbs free energy as

$$\mathcal{G}([\mathbf{Z}]) = \mathcal{G}([\mathbf{Z}]) - \mu_{Y_p} \cdot [\mathbf{m}], \qquad (4.112)$$

by using Eq. (4.52), and defining the concentration of the moieties as

$$[\mathbf{m}] = (\mathbb{1}_{Y_p}^b)^{-1}\mathbb{1}^b[\mathbf{Z}]. \qquad (4.113)$$

Equation (4.112) is reminiscent of the grand canonical free energy in equilibrium thermodynamics: when passing from the canonical to the grand canonical ensemble, the grand canonical free energy is obtained from the Gibbs free energy by eliminating the energetic contribution due to the matter exchange with the environment. Here, something similar happens: we eliminate the energetic contribution due to the exchanged moieties, i.e., $\mu_{Y_p} \cdot [\mathbf{m}]$. This is the reason why $\mathcal{G}([\mathbf{Z}])$ is also called a semi-grand Gibbs free energy. Note that $\mathcal{G}([\mathbf{Z}])$ is a state function and therefore its time derivative vanishes at steady-state. Note also that if $\dot{w}_{\text{nc}} = 0$ and $\dot{w}_{\text{driv}} = 0$, Eq (4.110) simplifies to

$$\frac{d}{dt}\mathcal{G}([\mathbf{Z}]) = -T\dot{\mathcal{S}}(t) \leq 0. \qquad (4.114)$$

Namely, the semigrand free energy is dissipated (decreasing monotonously in time) and CRNs relax towards an equilibrium state since $\mathcal{G}([\mathbf{Z}])$ is lower bounded by its equilibrium value $\mathcal{G}([\mathbf{Z}]_{\text{eq}})$ (see 4.6.6).

The nonconservative work rate is given by

$$\dot{w}_{\text{nc}} = \mathcal{F}_{Y_f} \cdot \mathbf{I}^{Y_f}, \qquad (4.115)$$

and quantifies the energetic cost to sustain fluxes of moieties between different chemostats through the CRN by means of the fundamental nonconservative forces

$$\mathcal{F}_{Y_f} = (\mu_{Y_f} \cdot \mathbb{1} - \mu_{Y_p} \cdot (\mathbb{1}_{Y_p}^b)^{-1}\mathbb{1}_{Y_f}^b)^T \qquad (4.116)$$

Indeed, $\dot{w}_{\text{nc}}$ always vanishes when all the chemostatted species break a conservation law, i.e., $S_Y = S_{Y_p}$. Only when there are some $S_{Y_f}$ species (i.e., when at least one moiety is exchanged with more than one chemostat), the forces (4.116) breaking the detailed balance condition and keeping the system out of equilibrium can emerge. This is the reason why the $\{Y_\alpha\}$ species with $\alpha \in S_{Y_f}$ are named force species, while the $\{Y_\alpha\}$ species with $\alpha \in S_{Y_p}$ are named potential species. Note also that the forces (4.116) can vanish, and hence the CRN becomes detailed balanced, even if there are force species. Indeed, $\mathcal{F}_{Y_f} = 0$ when $\mu_{Y_f} \cdot \mathbb{1} = \mu_{Y_p} \cdot (\mathbb{1}_{Y_p}^b)^{-1}\mathbb{1}_{Y_f}^b$. Thus, the fundamental nonconservative forces are generated by the difference between the actual values of chemical potentials of the force species,
i.e., $\mu_{Y_f}$, and the equilibrium values they would have if only the potential species were chemostatted, i.e., $\mu_{Y_p} \cdot (\mathbb{L}_{Y_p}^b)^{-1}\mathbb{L}_{Y_f}^b$.

*Remark.* The nonconservative forces $\mathcal{F}_{Y_f}$ and corresponding currents $\mathbf{I}^{Y_f}$ do not have to be aligned. This allows for free energy transduction.

The driving work rate $\dot{w}_{\text{driv}}$ is specified as

$$ \dot{w}_{\text{driv}} = -\frac{d}{dt}\mu_{Y_p} \cdot [\mathbf{m}], \qquad (4.119) $$

and accounts for the time dependent control of only the chemical potentials $\mu_{Y_p}$. It represents the energetic cost of modifying the equilibrium that is defined in the absence of nonconservative forces and vanishes in autonomous systems.

### 4.6.6 Thermodynamic potential

We conclude here by showing that $\mathcal{G}([\mathbf{Z}])$ is always lower bounded by its equilibrium value $\mathcal{G}([\mathbf{Z}]_{\text{eq}})$. To do so, we use a similar strategy to the one we used in Subs. (4.5.3).

First, we use that $\mu([\mathbf{Z}]_{\text{eq}})$ is still a linear combination of conservation laws (given in Eq. (4.72)) because of Eq. (4.44). However, only the unbroken conservation laws still define conserved quantities. Thus,

$$ \mu([\mathbf{Z}]_{\text{eq}}) \cdot [\mathbf{Z}]_{\text{eq}} = \sum_{\lambda_u} f_{\lambda_u} \underbrace{\ell^{\lambda_u} \cdot [\mathbf{Z}]_{\text{eq}}}_{=L^{\lambda_u}} + \sum_{\lambda_b} f_{\lambda_b} \ell^{\lambda_b} \cdot [\mathbf{Z}]_{\text{eq}} \qquad (4.120a) $$

$$ = \sum_{\lambda_u} f_{\lambda_u} \underbrace{\ell^{\lambda_u} \cdot [\mathbf{Z}]}_{=L^{\lambda_u}} + \sum_{\lambda_b} f_{\lambda_b} \ell^{\lambda_b} \cdot [\mathbf{Z}]_{\text{eq}} \qquad (4.120b) $$

$$ = \mu([\mathbf{Z}]_{\text{eq}}) \cdot [\mathbf{Z}] + \sum_{\lambda_b} f_{\lambda_b} \ell^{\lambda_b} \cdot [\mathbf{Z}]_{\text{eq}} - \sum_{\lambda_b} f_{\lambda_b} \ell^{\lambda_b} \cdot [\mathbf{Z}] \qquad (4.120c) $$

where we summed and subtracted $\sum_{\lambda_b} f_{\lambda_b} \ell^{\lambda_b} \cdot [\mathbf{Z}]$.

Second, we use that the entries of the unbroken conservation laws corresponding to the exchanged species vanish:

$$ \mu_{Y_p} = \sum_{\lambda_u} f_{\lambda_u} \underbrace{\ell_{Y_p}^{\lambda_u}}_{=0} + \sum_{\lambda_b} f_{\lambda_b} \ell_{Y_p}^{\lambda_b} = \mathbf{f}_b \cdot \mathbb{L}_{Y_p}^b, \qquad (4.121) $$

where $\mathbf{f}_b := (\dots, f_{\lambda_b}, \dots)^{\text{T}}$. This implies that

$$ \mu_{Y_p} \cdot [\mathbf{m}] = \mathbf{f}_b \cdot \mathbb{L}_{Y_p}^b [\mathbf{m}] \qquad (4.122) $$

$$ = \mathbf{f}_b \cdot \mathbb{L}_{Y_p}^b (\mathbb{L}_{Y_p}^b)^{-1} \mathbb{L}^b [\mathbf{Z}] \qquad (4.123) $$

$$ = \mathbf{f}_b \cdot \mathbb{L}^b [\mathbf{Z}] \qquad (4.124) $$

$$ = \sum_{\lambda_b} f_{\lambda_b} \ell^{\lambda_b} \cdot [\mathbf{Z}] \qquad (4.125) $$

$$ \mu_{Y_f}^{\text{eq}} = \sum_{\lambda_u} f_{\lambda_u} \underbrace{\ell_{Y_f}^{\lambda_u}}_{=0} + \sum_{\lambda_b} f_{\lambda_b} \ell_{Y_f}^{\lambda_b} = \mathbf{f}_b \cdot \mathbb{L}_{Y_f}^b, \qquad (4.117) $$

with $\mathbf{f}_b := (\dots, f_{\lambda_b}, \dots)^{\text{T}}$. This together with Eq. (4.121) and the existence of $(\mathbb{L}_{Y_p}^b)^{-1}$, namely, $\mu_{Y_p} \cdot (\mathbb{L}_{Y_p}^b)^{-1} = \mathbf{f}_b \cdot \mathbb{I}$ leads to

$$ \mu_{Y_f}^{\text{eq}} = \mu_{Y_p} \cdot (\mathbb{L}_{Y_p}^b)^{-1} \mathbb{L}_{Y_f}^b, \qquad (4.118) $$

as stated in the main text.

<footnote n="15">¹⁵To prove that $\mu_{Y_p} \cdot (\mathbb{L}_{Y_p}^b)^{-1}\mathbb{L}_{Y_f}^b$ is the equilibrium value of the chemical potentials of the force species if only the potential species were chemostatted, we use that the entries of the unbroken conservation laws corresponding to the exchanged species vanish and thus</footnote>
where we used the definition of the concentrations of the moieties given in Eq. (4.113), and analogously

$$ \mu_{Y_p} \cdot [m_{eq}] = \sum_{\lambda_b} f_{\lambda_b} \ell^{\lambda_b} \cdot [Z]_{eq}. \qquad (4.126) $$

Third, we use Eqs. (4.120c), (4.125) and (4.126) in the difference $G([Z]) - G([Z]_{eq})$ and we obtain

$$ G([Z]) - G([Z]_{eq}) = D([Z]) \| [Z]_{eq}) \ge 0, \qquad (4.127) $$

with $D([Z]) \| [Z]_{eq})$ defined in Eq. (4.76).

**Example. Autonomous and Detailed Balanced.** Consider the CRN (4.3) when only the species S and F are exchanged and their concentrations are maintained constants. Both species S and F are potential species. The semigrand free energy (4.112) thus reads

$$ G([Z]) = G([Z]) - \mu_F L^F - \mu_S L^S, \qquad (4.128) $$

where $\mu_F$ and $\mu_S$ are the (constant) chemical potential of F and S, respectively. For this setup, the second law is given in Eq. (4.114) since there are no force species and the concentrations [S] and [F] are constant. This, together with $G([Z])$ being lower bounded by its equilibrium value, implies that the CRN is detailed balance and will eventually reach an equilibrium state despite being open. Note that, by using the thermodynamic theory we derived, we can predict that the CRN will reach an equilibrium state without solving the dynamics. The only information we need is encoded in the stoichiometric matrix, i.e., in the topology of the CRN.

**Example. Nonautonomous and Detailed Balanced.** Consider the CRN (4.3) when only the species S and F are exchanged and their concentrations are changed in time according to unspecified protocols. The semigrand free energy (4.112) is still given in Eq. (4.128), since there are no new potential species, but the second law now becomes

$$ T \dot{\Sigma} = -\frac{d}{dt} G([Z]) + \dot{w}_{\text{driv}}, \qquad (4.129) $$

where

$$ \dot{w}_{\text{driv}} = -\left(\frac{d}{dt}\mu_F\right)L^F - \left(\frac{d}{dt}\mu_S\right)L^S. \qquad (4.130) $$

The CRN is still detailed balance, but the time dependent manipulation of the concentrations [F] and [S] provides free energy, in the form of the driving work (4.130), which changes the equilibrium to which the CRN would relax and balances the dissipation so maintaining the CRN out of equilibrium.

**Example. Autonomous and Nondetailed Balanced.** Consider the CRN (4.3) when the species P, W, S and F are exchanged and their concentrations are maintained constant. Exchanging the species P and W does not break any additional conservation laws, but creates the following nonconservative forces (obtained by specializing Eq. (4.116)):

$$ \mathcal{F}_P = \mu_P - \mu_S \quad \text{and} \quad \mathcal{F}_W = \mu_W - \mu_F. \qquad (4.131) $$

Hence, the semigrand free energy (4.112) is still given in Eq. (4.128), since there are no new potential species, but the second law now becomes

$$ T \dot{\Sigma} = -\frac{d}{dt} G([Z]) + \dot{w}_{\text{nc}}, \qquad (4.132) $$

with

$$ \dot{w}_{\text{nc}} = (\mu_P - \mu_S)I_P + (\mu_W - \mu_F)I_W. \qquad (4.133) $$

The CRN is not detailed balance because of the nonconservative forces (4.131) and will reach a nonequilibrium steady state. Note again that, by using the thermodynamic theory we derived, we can predict that the CRN will reach a nonequilibrium steady state without solving the dynamics.
4.7 Conclusions

We derived here, in a self-contained way, a thermodynamic theory for deterministic open CRNs in ideal dilute solutions. Crucially, we used the conservation laws to derive the formulation of the second law in Eq. (4.110) which identifies the specific mechanisms (nonconservative forces and driving) that can balance dissipation and maintain open CRNs out of equilibrium.

The same formulation of the second law (4.110) can be derived for deterministic non-
ideal CRNs [30], deterministic reaction-diffusion systems [69, 70], and stochastic CRNs in
ideal dilute solutions [29].
# 5 Stochastic Chemical Reaction Networks

**Matteo Polettini.** I illustrate results about the stochastic dynamics and thermodynamics of mass-action kinetics chemical reaction networks by simple examples and outputs of simulations. In particular I focus on the role of deficiency for stationary fluctuations.

## 5.1 Introduction

This lecture is an attempt to provide an incremental bottom-up presentation of techniques, results and nomenclature to treat the stochastic dynamics and thermodynamics of mass-action kinetics (MAK) chemical reaction networks (CRNs). The topics addressed would be enough for at least a full master course ranging from linear algebra to metabolic networks and quite possibly machine learning (ça va sans dire). In the two-hours span that was gently offered by the organizers I decided to provide by examples some basic intuitions, leaving generalizations and proofs to the willingness of the student and to the literature.

## 5.2 Fishing for Poisson

Consider the chemical reaction

$$ \emptyset \rightleftharpoons X. \qquad (5.1) $$

Let $\hat{X}(t)$ be the population (number of molecules present in the reactor) of species X at time $t^{16}$. Apart from species $\emptyset$ that has fixed population [$\emptyset](t) = 1$ (possibly because it is chemostatted via osmosis or other mechanisms, see Lecture 4), this is a fluctuating quantity taking integer values $X \in \mathbb{N}$, due to what is sometimes called intrinsic or population noise (to distinguish it from outside disturbances). We can represent the population space as

$$ 0 \text{ --- } 1 \text{ --- } 2 \text{ --- } 3 \text{ --- } \dots \qquad (5.2) $$

MAK prescribes that the rate of injection of molecules is constant, while that of the ejection of molecules is proportional to the population

$$ \begin{aligned} r_{+}(X) &= k_{+} \\ r_{-}(X) &= k_{-} X. \end{aligned} \qquad (5.3) $$

These are called reaction rates or velocities, while the $k_{\pm}$ are called rate constants.

Consider the probability $p_t(X) = p(\hat{X}(t) \equiv X)$ that at time $t$ the reactor contains $X$ molecules and define the (average) population-space current from $X-1$ to $X$ as

$$ j_t(X, X-1) = r_+(X-1)p_t(X-1) - r_-(X)p_t(X), \qquad (5.4) $$

with $j_t(0, -1) = 0$ The evolution of the probability is dictated by the Chemical Master Equation

$$ \frac{d}{dt} p_t(X) = j_t(X, X-1) - j_t(X+1, X) \qquad (5.5) $$

given some initial distribution $p_0$.

<footnote n="16">¹⁶In this lecture, differing from Lecture 4, we use capital italic letters X for stochastic populations and small-case italic letters x for deterministic concentrations in place of [X]).</footnote>
from random import uniform

from math import log

kp = 2 # rate constant forward

km = 1 # rate constant backward

x = 0 # initial population

t = 100000. # total time

tau = 0 # time elapsed

h1, h2 = {}, {} # output histograms

while tau < t:

rp = kp # calculate rate forward with MAK

rm = km * x # calculate rate backward with MAK

r = rp + rm

y1 = uniform(0,1)

y2 = uniform(0,1)

s = log(1/y1)/r # sample time by inversion rule

h1[x] = h1.get(x, 0) + s # count time spent at state

h2[x] = h2.get(x,0) + 1 # count state occurrences

x == x - 1 if y2 < rp / r else x == x + 1 # sample new state

tau += s

Figure 5.1: Doob-Gillespie algorithm to generate the histograms of the total number of occurrences of a population and of the total time spent at a state and compare to theory. The output is shown in Fig. 5.2.