Phase transitions in thermodynamically consistent biochemical systems Von der Fakultät Mathematik und Physik der Universität Stuttgart zur Erlangung der Würde eines Doktors der Naturwissenschaften (Dr. rer. nat.) genehmigte Abhandlung vorgelegt von Basile Nguyen aus Lausanne (Schweiz) Hauptberichter: Prof. Dr. Udo Seifert Mitberichter: Prof. Dr. Siegfried Dietrich Vorsitzender: Prof. Dr. Martin Dressel Tag der Einreichung: 05. Juni 2020 Tag der mündlichen Prüfung: 22. Juli 2020 II. Institut für Theoretische Physik der Universität Stuttgart 2020 Contents Publications 7 Zusammenfassung 9 Summary 15 1 Introduction 19 2 Basics 23 2.1 Stochastic thermodynamics for discrete state spaces . . . . . . . . . . . 23 2.2 Chemical reaction networks . . . . . . . . . . . . . . . . . . . . . . . . 27 3 First-order phase transitions in biochemical switches 41 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2 Schlögl’s model and entropy production . . . . . . . . . . . . . . . . . 42 3.3 Chemical master equation . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.4 Two-state model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.5 Fokker-Planck equation . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Appendices to chapter 3 61 3.A Relation between the cumulants of currents in a system with two reaction channels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.B Calculation of the coarse-grained rates for the two-state model . . . . . 61 3.C Calculation of the diffusion coefficient without relying on large deviation theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.D Tilted operator for the Fokker-Planck equation . . . . . . . . . . . . . . 69 4 Second-order phase transitions in biochemical oscillators 71 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.2 Brusselator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.3 Activator-inhibitor model . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.4 KaiC model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.5 Criteria for the precision of biochemical oscillations . . . . . . . . . . . 85 4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 3 Contents Appendices to chapter 4 91 4.A Activator-inhibitor model . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.B KaiC model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5 Concluding perspective 95 Bibliography 99 Danksagung 111 4 List of Figures 3.1 Phase transition in the Schlögl model . . . . . . . . . . . . . . . . . . . 50 3.2 Diffusion coefficient for the Schlögl model . . . . . . . . . . . . . . . . 51 4.1 Trajectories at different chemical forces . . . . . . . . . . . . . . . . . 74 4.2 Phase transition in the Brusselator . . . . . . . . . . . . . . . . . . . . 75 4.3 Diffusion coefficient for the Brusselator . . . . . . . . . . . . . . . . . 76 4.4 Activator-inhibitor model . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.5 Phase transition for the activator-inhibitor model . . . . . . . . . . . . . 80 4.6 Diffusion coefficient for the activator-inhibitor model . . . . . . . . . . 81 4.7 KaiC model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.8 Phase transition for the KaiC model . . . . . . . . . . . . . . . . . . . 83 4.9 Diffusion coefficient for the KaiC model . . . . . . . . . . . . . . . . . 84 4.10 Correlation function for the oscillating species . . . . . . . . . . . . . . 86 4.11 Comparison between the number of coherent oscillations and the Fano factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5 Publications Parts of this thesis have been published in • B. Nguyen and U. Seifert, Exponential volume dependence of entropy-current fluctuations at first-order phase transitions in chemical reaction networks, Phys. Rev. E 102, 022101 (2020). © 2020 American Physical Society. Chapter 3 is based on this publication. • B. Nguyen, U. Seifert and A.C. Barato, Phase transition in thermodynamically consistent biochemical oscillators, J. Chem. Phys. 149, 045101 (2018). © 2018 AIP Publishing. Chapter 4 is based on this publication. Further publications I have co-authored • J. H. Fritz, B. Nguyen and U. Seifert, Stochastic thermodynamics of chemical reactions coupled to finite reservoirs: A case study for the Brusselator, J. Chem. Phys. 152, 235101 (2020). • A. Ehrmann, B. Nguyen, U. Seifert, Interlinked GTPase cascades provide a motif for both robust switches and oscilla- tors, J. R. Soc. Interface 16, 20190198 (2019). • B. Nguyen, D. Hartich, U. Seifert and P. De Los Rios, Thermodynamic bounds on the ultra- and infra-affinity of Hsp70 for its substrates, Biophys. J. 113, 362-370 (2017). 7 https://10.1103/PhysRevE.102.022101 https://doi.org/10.1063/1.5032104 https://doi.org/10.1063/5.0006115 https://doi.org/10.1098/rsif.2019.0198 https://doi.org/10.1016/j.bpj.2017.06.010 Zusammenfassung Lebende Systeme benötigen zur Erfüllung ihrer Aufgaben eine konstante Energiezufuhr. Viele biologische Systeme benötigen Energieträger, die sogenannten Nukleotidphos- phate, wie z.B. Adenosintriphosphat (ATP). Chemische Energie wird bei einer Hydro- lysereaktion durch Aufbrechen einer der Phosphatbindungen erzeugt. Diese Reaktion verwandelt ein ATP in ein Adenosindiphosphat (ADP) und ein Phosphat (Pi). Die Zellen halten einen ATP-Exzess durch einen metabolischen Prozess, die sogenannte Zellat- mung, aufrecht, wobei ADP kontinuierlich in ATP umwandelt wird. Zur Vereinfachung können wir davon ausgehen, dass ein biologisches System an chemische Reservoirs ge- koppelt ist, die viel größer sind als das System und dass sich diese Reservoirs auf den relevanten Zeitskalen daher nicht signifikant ändern. Wenn das System nicht von außen gesteuert wird, kommt es ab einem gewissen Zeitpunkt zu einem stationären Nicht- Gleichgewichtszustand. Biologische Systeme beruhen auf Bistabilität und Oszillationen, um entscheidende Prozesse zu regulieren. Einerseits steuern Schalter die zelluläre Reaktion auf ein Input- Signal über komplexe Rückkopplungsstrukturen und sind zentrale Elemente vieler zellu- lärer Prozessewie derGenexpression, desMembrantransports oder desVesikeltransports. Ein idealer Schalter sollte ein all-or-none Response aufweisen und auf Bistabilität beru- hen, um gegenüber Input-Störungen robust zu sein. Andererseits steuern biochemische Oszillatoren das Timing und die Kontrolle der zentralen Prozesse der Zelle wie zirkadia- ne Rhythmen und den Zellzyklus. Diese zwei Regime sind mit Phasenübergängen erster und zweiter Ordnung verbunden. In dieser Dissertation werden wir den Phasenübergänge charakterisieren, der in biochemischen Systemen stattfindet. Laut der Thermodynamik können komplexe Verhaltensweisen in chemischen Reakti- onsnetzwerken wie Bistabilität oder Oszillationen nur weit außerhalb des Gleichgewichts erreicht werden. ImGleichgewicht wird das System in einen eindeutigen Gleichgewichts- zustand relaxieren. Das Verhalten des Systems hängt von seinem Abstand vom Gleich- gewicht ab und muss an einem gewissen Punkt zu einem Phasenübergang kommen, um Oszillationen auszubilden. Bei Systemen, die in Kontakt mit chemischen Reservoirs sind, kann die Differenz der chemischen Potentiale zwischen den Reservoirs ein Steuer- parameter für diesen Phasenübergang sein. Diese Differenz führt zu thermodynamischen Kräften, die das System in einem Nicht-Gleichgewichtszustand treiben. Zum Beispiel hängt in Zellen die thermodynamische Kraft von der Differenz der chemischen Potentia- le zwischen ATP, ADP und Pi ab. 9 Zusammenfassung Nichtgleichgewichts-Phasenübergänge wurden seit langem untersucht und sind immer noch weniger verstanden als ihre Pendants im Gleichgewicht. Das Konzept der Ord- nungsparameter kann auf ähnliche Weise wie bei Gleichgewichtssystemen verwendet werden. Allerdings kann es nicht den irreversiblen Charakter eines Nichtgleichgewichts- Phasenübergangs beschreiben. Die neuesten Fortschritte in der stochastischen Thermo- dynamik haben die Wichtigkeit der Entropieproduktion als Quantifikator der Entfernung vom Gleichgewicht betont, da sie in einem Gleichgewichtssystem verschwindet und in einem Nicht-Gleichgewichtsregime positiv ist. Daraus folgt die Frage, ob die Entropie- produktion ein charakteristisches Verhalten bei Phasenübergängen zeigt. Die Erforschung von Stromfluktuationen hat sich seit der Entdeckung der “Thermody- namischen Unschärferelation” im Jahr 2015 wesentlich weiterentwickelt. Diese Relation besagt, dass die Präzision jedes Stroms einen minimalen thermodynamischen Kosten- faktor hat. Sie gilt für eine breite Palette von Systemen, die durch einen Markov-Prozess beschrieben werden. Ein besonders wichtiger Strom ist die totale Entropieprodukti- on. Obwohl die Entropie-Stromfluktuationen eine universelle untere Schranke haben, sind sie von oben nicht beschränkt, insbesondere in der Nähe eines Phasenübergangs, bei dem wir erwarten, dass die Fluktuationen thermodynamischer Größen divergieren. Das Hauptziel dieser Dissertation ist eine Charakterisierung von Nichtgleichgewichts- Phasenübergängen. Insbesondere interessieren wir uns für das Verhalten der Entropie- produktion bei Phasenübergängen erster und zweiter Ordnung. Kapitel 2: Grundlagen Dieses Kapitel bietet einen Überblick über die Grundlagen der stochastischen Dynamik, der stochastischen Thermodynamik und der chemischen Reaktionsnetzwerke. Ausgehend von einer grundlegenden Theorie der statistischen Mechanik beginnen wir mit einem ge- schlossenen System in Kontakt mit einem Wärmebad und identifizieren Mesozustände. Angenommen das System äquilibriert lokal innerhalb jedes Mesozustände, dann wird die Dynamik durch einen Markov-Prozess auf der mesoskopischen Ebene beschrieben. Ein geschlossenes System wird ein thermisches Gleichgewicht erreichen, daher müssen die Übergangsraten zwischen den Mesozuständen eine detaillierte Gleichgewichtsbedin- gung erfüllen, die die freien Energien der Mesozustände beinhaltet. Durch die Aufteilung eines geschlossenen Systems in ein Kernsystem und chemische Reservoirs können wir die Dynamik offener Systeme beschreiben. Durch entsprechende Steuerung der che- mischen Potenziale der Reservoirs kann ein offenes System im Gegensatz zu seinem geschlossenen Pendant einen stationären Nicht-Gleichgewichtszustand erreichen. Dann identifizieren wir thermodynamische Observablen wie Wärme-, Arbeits- oder Entropie- produktion entlang einzelner Trajektorien so, dass sich die erwarteten makroskopischen Observablen ergeben. Daraufhin betrachten wir die Dynamik chemischer Reaktionsnetzwerke. Mit der An- 10 nahme, dass chemische Reaktionen Poisson-Prozesse sind, definieren wir Übergangs- raten und führen die chemische Master-Gleichung (CME) ein, die die Dynamik eines chemischen Reaktionsnetzwerks beschreibt. Es ist im Allgemeinen sehr schwierig diese Master-Gleichung exakt zu lösen. Wir stellen daher den stochastischen Simulationsalgo- rithmus vor, der eine exakte Methode zur Erzeugung stochastischer Trajektorien liefert. Zunächst leiten wir die chemische Langevin-Gleichung (“CLE”) aus der CME ab, in- dem wir die Poisson-Zufallsvariablen mit normalverteilten Zufallsvariablen annähern. Die CLE kann auch als Fokker-Planck-Gleichung formuliert werden, die man normaler- weise durch Trunkierung der chemischen Kramers-Moyal-Gleichung erhält. Wir leiten die lineare Rauschnäherung nach van Kampen ab, indem wir die CLE linearisieren. Im Limes unendlicher Systemgröße erhalten wir die deterministischen Ratengleichungen. Schließlich diskutieren wir, ob trimolekulare Reaktionen sinnvoll sind. Wir zeigen, dass eine trimolekulare Reaktion auf der stochastischen Ebene die gleiche Dynamik haben kann wie eine Reihe bimolekularer Reaktionen, bei denen eine zusätzliche kurzlebige Spezies eingeführt werden muss. Kapitel 3: Phasenübergang erster Ordnung bei biochemischen Schaltern In diesem Kapitel betrachten wir biochemische Schalter, die bei ihrer Aktivierung einen Phasenübergang erster Ordnung durchlaufen. Aus einer deterministischen Perspektive ist dieser Phasenübergang mit dem Konzept der Bistabilität verbunden, bei der zwei stabile stationäre Zustände koexistieren können. Im Gegensatz dazu ist in der stochastischen Per- spektive der stationäre Zustand eindeutig und mit einer bimodalen Verteilung verbunden. Bei Phasenübergängen erster Ordnung ist die Entropieproduktionsrate bekanntermaßen diskontinuierlich. Wir zeigen, dass die Fluktuationen der Entropieproduktion in der Nähe des bistabilen Punktes eine exponentielle Volumenabhängigkeit haben. Unsere Ergeb- nisse, die wir für ein paradigmatisches Modell biochemischer Schalter, das Schlögl’sche Modell, erhalten haben, können auf eine breite Klasse von Nichtgleichgewichtssystemen angewendet werden. Wir beginnen mit der chemischen Master-Gleichung. Mit der Theorie der großen Abweichungen berechnen wir die mittlere Rate der Entropieproduktion und den zuge- hörigen Diffusionskoeffizienten numerisch. Im Limes großer Systemgrößen gibt es zwei relevante Zeitskalen: eine schnelle Relaxation zum nächsten Fixpunkt und einen langsa- meren Übergang zwischen den Fixpunkten. In der Nähe der Fixpunkte kann das System durch stationäre Gauß’sche Prozesse modelliert werden. Für die langsamere Dynamik reduzieren wir die CME auf ein Zwei-Zustand Modell, in dem wir zeigen, dass die ver- gröberten Übergangsraten proportional zum Exponential des inversen Volumens sind. Für dieses Zwei-Zustands-Modell berechnen wir einen analytischen Ausdruck für den Diffusionskoeffizienten und zeigen, dass er eine exponentielle Volumenabhängigkeit am 11 Zusammenfassung bistabilen Punkt besitzt, wobei der exponentielle Vorfaktor durch die Höhe der effektiven Potentialbarriere, die die beiden Fixpunkte trennt, gegeben ist. Die Berechnung der Diffusionskoeffizienten für ein Zwei-Zustands-Modell ist relativ einfach, die Lösung wird jedoch in Form der vergröberten Übergangsraten ausgedrückt. Die Ableitung dieser Raten beginnend mit der CME kann eine anspruchsvolle Aufgabe sein, zu sehen im Appendix dieses Kapitels. In einem alternativen Zugang beginnend mit der Fokker-Planck-Gleichung leiten wir einen analytischen Ausdruck für den Diffusi- onskoeffizienten her. Für bistabile Systeme beschreibt die zweifach verkürzte chemische Kramers-Moyal-Gleichung nur die Dynamik auf der kürzeren Zeitskala richtig und hat außerdem nicht die richtige stationäre Verteilung. Hier betrachten wir die Fokker-Planck- Gleichung mit einem effektiven Diffusionsterm, der zu derselben stationären Verteilung wie die CME führt und die Dynamik auf der längeren Zeitskala korrekt beschreibt. Kapitel 4: Phasenübergang zweiter Ordnung in biochemischen Oszillatoren Biochemische Oszillationen sind in lebenden Organismen omnipräsent. In einem au- tonomen System, das nicht durch ein externes Signal beeinflusst wird, können sie nur im Nichtgleichgewicht auftreten. Wir zeigen, dass sie durch einen generischen Nicht- Gleichgewichts-Phasenübergang entstehen, mit einem charakteristischen qualitativen Verhalten bei Kritikalität. Wir untersuchen das kritische Verhalten des Entropieprodukti- onsstroms für drei Versionen bekannter Modelle: das Brusselator-Modell, das Aktivator- Inhibitor-Modell und ein Modell für KaiC-Oszillationen. Wir stellen fest, dass die erste Ableitung der mittleren Rate der Entropieproduktion eine Diskontinuität am kritischen Punkt besitzt, an dem die Fluktuationen als Potenzgesetz mit dem Volumen divergieren. Anschließend werden Metriken für die Präzision der biochemischen Oszillationen diskutiert. Als erstes betrachten wir das Standardmaß für die Präzision, die Anzahl der kohärenten Oszillationen, d.h. die Anzahl der Perioden, für die verschiedene stochasti- sche Realisationen miteinander kohärent sind. Zweitens betrachten wir den Fano-Faktor, der mit dem thermodynamischen Fluss verbunden ist. Angesichts der thermodynami- schen Unschärferelation ist der Fano-Faktor eine logische Wahl, um die Präzision mit der Thermodynamik in Beziehung zu setzen. Außerdem wurde er als Messwert vorge- schlagen, um die Präzision der biochemischen Oszillatoren zu quantifizieren. Da der Fano-Faktor selbst dann klein sein kann, wenn es keine biochemischen Oszillationen gibt, argumentieren wir, dass die Anzahl der kohärenten Oszillationen besser geeignet ist, um die Präzision der biochemischen Oszillationen zu quantifizieren. 12 Ausblick Das Hauptziel dieser Arbeit war es, das Verhalten des Entropieproduktionsstroms und seiner Fluktuationen bei Phasenübergängen erster und zweiter Ordnung in homogenen chemischen Reaktionsnetzwerken zu beschreiben. Für chemische Reaktionsnetzwerke mit räumlicher Dynamik, in denen chemische Reaktionen lokal stattfinden und sich die chemische Spezies durch Diffusion ausbreitet, bleibt die Charakterisierung der Entropie- produktion eine offene Frage. Überraschenderweise bleibt die gesamte Entropieproduk- tion in solchen Systemen beim Phasenübergang analytisch. Wir haben die Präzision der biochemischen Oszillationen in autonomen Systemen untersucht. Eine interessante zukünftige Forschungsrichtung wären biochemische Os- zillatoren, die von einem externen periodischen Signal angetrieben werden. Für solche Systeme wurde gezeigt, dass die Präzision der Oszillationen unendlich sein kann. Au- ßerdem hat sich in letzter Zeit gezeigt, dass Stromfluktuationen verschwinden können und damit die ursprüngliche thermodynamische Unschärferelation verletzt wird. Eine weitere unerwartete Möglichkeit, die Präzision der Oszillationen zu verbessern, ist die Berücksichtigung endlicher chemischer Reservoirs. Wir haben kürzlich gezeigt, dass sie ihre idealen Pendants übertreffen können. 13 Summary Living systems require a constant supply of energy to perform their tasks and they must operate far from equilibrium. Most biological systems rely on energy carriers called nucleotide phosphates such as adenosine triphosphate (ATP). Chemical energy can be extracted through a hydrolysis reaction by breaking one of the phosphate bonds. This reaction converts an ATP into an adenosine diphosphate (ADP) and an inorganic phosphate (Pi). Cells maintain an excess of ATP through a metabolic process called cellular respiration that continuously recycles ADP into ATP. For simplicity, we can assume that a biological system is coupled to chemical reservoirs that are much larger than the system and thus these reservoirs do not evolve in a significant way on relevant timescales. If the system is not externally controlled, it will at some point in time reach a stationary state, which is called a nonequilibrium steady-state (NESS). Biological systems rely on bistability and oscillations to regulate key processes. On one hand, switches control the cellular response to an input signal using complex feed- back structures and are key players in many cellular processes such as gene expression, membrane trafficking or vesicular transport. An ideal switch should have an all-or-none response and rely on bistability to be robust against input noise. On the other hand, biochemical oscillators set the timing and control of key processes such as circadian rhythms and the cell cycle. These two regimes are associated with first- and second-order phase transitions, respectively. In this thesis, we will characterize the phase transition occurring in biochemical systems. Thermodynamics states that complex behaviors in chemical reactions networks such as bistability or oscillations can only be achieved far from equilibrium. In equilibrium, the system will relax towards a unique equilibrium steady-state. The behavior of the system depends on its distance from equilibrium, and thus it must undergo a phase transition at some point. For systems in contact with chemical reservoirs, a control parameter for this phase transition can be a difference of chemical potentials between the reservoirs resulting in thermodynamic forces driving the system in a NESS. For instance, in cells the thermodynamic force depends on the difference in chemical potentials between ATP, ADP and Pi. Nonequilibrium phase transitions have long been studied and still remain less under- stood than their equilibrium counterparts. The concept of order parameter can be used in a similar way to equilibrium systems. However, it will not describe the irreversible character of a nonequilibrium phase transition. Recent advances in stochastic thermo- 15 Summary dynamics have emphasized the importance of entropy production as a quantifier for the distance from equilibrium as it vanishes for equilibrium systems and is positive in a nonequilibrium regime. A natural question that arises is whether the entropy production displays a characteristic behavior at phase transitions. The study of current fluctuations has recently seen some major developments since the discovery of the “thermodynamic uncertainty relation” in 2015. This relation states that the precision of any current has a minimal thermodynamic cost and holds for a large class of systems described by a Markov process. A current of particular importance is the total entropy production. While entropy current fluctuations have a universal lower bound, they are not bounded from above, especially in the vicinity of a phase transition where we expect the fluctuations of thermodynamic quantities to diverge. The main objective of the work outlined in this thesis is to provide a comprehensive characterization of nonequilibrium phase transitions. In particular, we are interested in the behavior of the entropy production at first- and second-order phase transitions. Chapter 2: Basics We review the basic principles of stochastic dynamics, stochastic thermodynamics and chemical reaction networks. Assuming an underlying theory of statistical mechanics, we start with a closed system in contact with a heat bath and identify mesostates. As- suming that the system equilibrates locally within each mesostate, the dynamics become Markovian at the mesoscopic level. A closed system will reach thermal equilibrium, thus, the transition rates between the mesostates must fulfill a detailed balance condition that involves the free energies of the mesostate. By separating the closed system into a core system and chemical reservoirs, we can describe the dynamics of open systems. By tuning the chemical potentials of the reservoirs accordingly, an open system can achieve a NESS contrary to its closed counterpart. We then identify thermodynamic observables such as heat, work or entropy production at the level of individual trajectories in such a way that the expected macroscopic observables are recovered. We then review the dynamics of chemical reaction networks. With the assumption that chemical reactions are Poisson processes, we define transition rates and introduce the chemical master equation (CME), which describes the dynamics of a chemical reaction network. In general, it is very difficult to solve this master equation exactly. We thus introduce the stochastic simulation algorithm, an exact procedure to generate stochastic trajectories. We then present a series of approximations of the CME. First, we derive the chemical Langevin equation (CLE) from the CME by approximating Poisson random variables with normal random variables. The CLE can also be expressed as Fokker- Planck equation usually obtained by truncating the chemical Kramers-Moyal equation. We derive the van Kampen’s linear noise approximation by linearising the CLE. In the limit of infinitely large system sizes, we obtain the deterministic rate equations. 16 Finally, we discuss whether trimolecular reactions are physically justified. We show that a trimolecular reaction can have the same dynamics at the stochastic level as a set of bimolecular reactions, where an additional short-lived species must be introduced. Chapter 3: First-order phase transition in biochemical switches In this chapter, we focus on biochemical switches that undergo a first-order phase transi- tion upon activation. From a deterministic perspective, this phase transition is associated with bistability where two stable steady states can coexist. In contrast, in the stochastic perspective, the steady-state is unique and is associated with a bimodal density dis- tribution. At first-order phase transitions, the entropy production rate is known to be discontinuous. We show that the fluctuations of the entropy production have an exponen- tial volume-dependence in the vicinity of the bistable point. Our results obtained for a paradigmatic model of biochemical switches, Schlögl’s model, can be applied to a large class of nonequilibrium systems. We start with the chemical master equation. Using large deviation theory, we compute the mean rate of entropy production and its associated diffusion coefficient numerically. In the limit of large system sizes, there are two relevant timescales: a fast relaxation to the nearest fixed-point and a slower transition between the fixed-points. Close to the fixed-points, the system can be modeled by stationary Gaussian processes. For the slower dynamics, we coarse-grain the CME into a two-statemodel, wherewe show that transition rates are proportional to the exponential of the inverse volume. For this two-state model, we compute an analytical expression for the diffusion coefficient and show that it has an exponential volume-dependence at the bistable point, where the exponential prefactor is given by the height of the effective potential barrier separating the two fixed-points. Computing diffusion coefficients for a two-state model is straightforward, however, the solution is expressed in terms of the coarse-grained transition rates. The derivation of these rates starting with CME can be challenging, as shown in the appendices of this chapter. As an alternative approach, we derive an analytical expression for the diffusion coefficient starting with the Fokker-Planck equation. For bistable systems, the two-term truncated chemical Kramers-Moyal equation does only correctly describe the dynamics on the shorter timescale and, moreover, does not have the correct stationary distribution. In our case, we consider a Fokker-Planck equation with an effective diffusion term that leads to the same stationary distribution as the CME and correctly describes the dynamics on the larger timescale. 17 Summary Chapter 4: Second-order phase transition in biochemical oscillators Biochemical oscillations are ubiquitous in living organisms. In an autonomous system, not influenced by an external signal, they can only occur out of equilibrium. We show that they emerge through a generic nonequilibrium phase transition, with a characteristic qualitative behavior at criticality. We characterize the critical behavior of the entropy production current for three versions of known models: the Brusselator, the activator- inhibitor model, and a model for KaiC oscillations. We find that the first derivative of the mean rate of entropy production has a discontinuity at the critical point, where the fluctuations diverge as a power-law with the volume. We then discuss metrics for the precision of biochemical oscillations. First, we consider the standard measure of precision, the number of coherent oscillations, which is the number of periods for which different stochastic realizations remain coherent with each other. Second, we consider the Fano factor associated with the thermodynamic flux. In light of the thermodynamic uncertainty relation, the Fano factor is a natural observable to relate precision with thermodynamics. It has been proposed as an observable that can quantify the precision of biochemical oscillators. Since the Fano factor can be small even when there are no biochemical oscillations, we argue that the number of coherent oscillations is more appropriate to quantify the precision of biochemical oscillations. Concluding perspective The main goal of this thesis was to describe the behavior of the entropy production current and its fluctuations at first- and second-order phase transitions in homogenous chemical reactions networks. For chemical reaction networks with spatial dynamics, where chemical reactions occur locally and diffusion spreads the chemical species, the characterization of the entropy production remains an open question. Surprisingly, the total entropy production in such systems remains analytical at the phase transition. We investigated the precision of biochemical oscillations in autonomous systems. As a future direction, it would be interesting to study biochemical oscillators under the influence of an external periodic signal. For these systems, the precision of oscillations can, in principle, be infinite. Moreover, it has been shown that current fluctuations can vanish, thus, violating the original thermodynamic uncertainty relation. Another unexpected way to improve the precision of oscillations is to consider finite reservoirs. We have recently found that they can outperform their ideal counterparts. 18 1 Introduction Thermodynamics is a universal theory for describing exchange processes and started with the development of heat engines that can deliver mechanical work using heat reservoirs. In 1824, Carnot proposed a theoretical engine consisting of a set of ideal operations, the Carnot cycle. This engine transfers heat from a warm region to a cold region, converting some of this heat into mechanical work [1]. Thanks to experiments by Joule, the equivalence of heat andmechanical energywas established, and around 1850Clausius and Thomson stated the two laws of thermodynamics. The first law is the principle of the conservation of energy, it relates the change in internal energy of a system with the extracted work and heat exchanged with the environment. The second law as formulated by Clausius states that “heat can never pass from a colder to a warmer body without some other change, connected therewith, occurring at the same time”. In an irreversible cycle, the heat transferred is always greater than the extracted work. This uncompensated heat is now called the total change in entropy production. Standard thermodynamics is limited to macroscopic systems in equilibrium. Onsager established a first formulation of nonequilibrium thermodynamics in the linear regime, i.e for small perturbations close to equilibrium. He derived universal symmetries, his recip- rocal relations, which couples currents to thermodynamic forces with phenomenological coefficients [2, 3]. The connection between thermodynamics and chemical reaction net- work was pioneered by Prigogine who proposed the concept of local equilibrium that enables the description of irreversible processes using equilibrium quantities [4]. In the last twenty years, nonequilibrium thermodynamics has substantially evolved with stochastic thermodynamics [5–7]. This framework extends the concepts of standard ther- modynamics such as work, heat or entropy to the level of individual trajectories. It can describe a large class of small fluctuating systems coupled to a heat bath and arbitrarily far from equilibrium. One of the major successes of this theory is the description of bi- ological processes. For instance, it has been successfully applied to cytoskeletal motors such as kinesin or myosin [8–11] or the F1-ATPase rotary motor [12, 13]. In the latter case, experimental techniques have enabled the observations of individual trajectories at the single-molecule level [14–16]. From such trajectories, efficiencies for the F1-ATPase could be computed [17–19]. Stochastic thermodynamics has recently seen major developments in understanding the cost of precision in nonequilibrium steady states. At its origin the discovery of the thermodynamic uncertainty relation (TUR) by Barato and Seifert [20] in 2015, it 19 1 Introduction was later proved by Gingrich and Horowitz [21] in 2017. This relation is an inevitable trade-off between the cost and precision of any current in a thermodynamically consistent system. In other words, small fluctuations can only be achieved by dissipating a minimal amount of entropy production. From a practical perspective, the TUR can be used to infer hidden properties of complex biochemical networks by analyzing experimental data from single-molecule experiments, in particular, the total entropy production. Most remarkably, the TUR implies an upper bound on the efficiencies of molecular motors [22, 23]. Furthermore, it states that a steady-state heat engine reaching Carnot efficiency at finite power must have diverging power fluctuations [24–27]. From a theoretical standpoint, one possibility to realize such an efficient engine would be to operate it in the vicinity of a phase transition, where fluctuations of observables are known to diverge. In the work outlined in this thesis, we will provide a comprehensive characterization of nonequilibrium phase transitions. We are particularly interested in the behavior of the entropy production current and its fluctuations at first- and second-order phase transitions. Phase transitions are arguably one of the most interesting phenomena in the field of physics, they have been long studied in the field of thermodynamics since its origin and remain an active field of research in modern statistical physics. Equilibrium phase transitions have a rich history and arewell understood from a thermodynamical standpoint [28, 29]. The first classification scheme for the transition between different phases of matter was introduced by Paul Ehrenfest [30] in 1933 after Keesom and Clusiys [31] discovered a phase transition in liquid helium in 1932, the lambda transition. Ehrenfest proposed that the order of a phase transition to be given by the lowest derivative of the free energy, which is discontinuous at the transition. For instance, the liquid-gas transition would classify as a first-order phase transition as it has a discontinuous change in density, which is the inverse of the first derivative of the free energy with respect to pressure. However, Onsager’s solution of the two-dimensional Ising model [32] had a free energy derivative with a power-law divergence rather than a discontinuity close to the critical point. This was experimentally confirmed in experiments by Atkins and Edwards who showed that the lambda transition was of a similar type [33]. This phase transition, which inspired Ehrenfest, turned out to be outside of his classification scheme. By the seventies, the classification scheme was simplified to binary classification with first-order and continuous transitions [34]. Nonequilibrium phase transitions have been long studied since the seventies [35–40] but still remain less understood than their equilibrium counterparts from a thermody- namical perspective. They have been originally classified using nonlinear dynamics and take place at bifurcations [41]. These bifurcations occur in biological systems that can achieve rich dynamics such as biochemical switching and oscillations, which are asso- ciated with saddle-node and Hopf bifurcation, respectively. These complex behaviors can be observed, for instance, in interlinked GTPases [42–45] or in the MinDE system [46–50], these generic structures can be tuned be as switches or oscillators depending on 20 their tasks. Their complex behavior can be understood with simple chemical networks introduced by Schlögl in the seventies to study nonequilibrium first- and second-order phase transitions [51]. Subsequently, order parameters with their associated variances have been derived for these systems [35, 39, 52–55]. Concepts from equilibrium systems such as the order parameter can be used to under- stand nonequilibrium phase transitions. However, these quantities cannot describe their irreversible character. For instance, the entropy production rate is a quantifier for the distance from equilibrium. It vanishes in equilibrium and is positive in a nonequilibrium regime [56]. A natural question that arises is whether the entropy production rate displays a characteristic behavior at phase transitions. This will be the main focus of this work. In this thesis, we study first- and second-order phase transitions in chemical reaction networks. We especially focus on the behavior of the entropy production current and its fluctuations. In Chapter 2, we review the basic principles of stochastic dynamics, stochastic thermodynamics, and chemical reaction networks. In Chapter 3, we investigate a paradigmatic model of biochemical switches, Schlögl’s model, that undergoes a first- order phase transition upon activation. We start with the chemical master equation and show numerically that the diffusion coefficient diverges at the bistable point. By coarse- graining the chemical master equation, we obtain a two-state model, where we show that the transition rates are proportional to the inverse of the system size. We obtain an analytical expression for the diffusion coefficient and show that it has an exponential volume dependence at the critical point, where the exponential prefactor is the height of the effective potential barrier separating the two fixed-points. As an alternative approach, we also compute the diffusion coefficient starting with the Fokker-Planck equation. In Chapter 4, we study a generic second-order phase transition occurring in autonomous biochemical oscillators. We characterize the critical behavior of the entropy production current in three known models: the Brusselator, the activator-inhibitor model and a model for KaiC oscillations. We show that the first derivative of the mean rate of entropy production is discontinuous at the critical point and that fluctuations diverge as a power- law with the volume. We discuss metrics for the precision of biochemical oscillations by comparing two observables, the Fano factor associated with the thermodynamic flux and the number of coherent oscillations. Since the Fano factor can be small in the absence of oscillations, we argue that the number of coherent oscillations is a better quantifier. In addition, we identify a region where the number of coherent oscillations and the current fluctuations both increase. Finally, we conclude in Chapter 5. 21 2 Basics 2.1 Stochastic thermodynamics for discrete state spaces In this section, we review the concepts of stochastic thermodynamics in discrete state space following [7, 57]. 2.1.1 Closed systems in contact with a heat bath We start with a closed system weakly coupled to heat bath at inverse temperature β. The system is fully described by its microstate ξ that contains the positions and momenta of all N particles in the system. Each microstate has an energy H(ξ), which may not be known explicitly. In equilibrium, the probability to observe a microstate is given by pcan(ξ) = exp [−βH(ξ)] Z(β) , (2.1) where the normalization constant is called the canonical partition function Z(β) ≡ ∫ dξ h3N exp [−βH(ξ)] . (2.2) Here, h is an an arbitrary quantity with the dimension of an action. The free energy, internal energy and entropy of the system are given by F = −(1/β) lnZ(β), U = ∂ ∂β (βF) , S(β) = β2 ∂ ∂β F = β (U − F) , (2.3) respectively. Solving the full equation of motion at the microscopic scale is neither possible nor of interest. We group the microstates ξ into a mesostate n and describe the dynamics at the mesoscopic level. Furthermore, we assume that each microstate ξ belongs to exactly one mesostate n. In equilibrium, the probability to observe a mesostate n is Peq(n) = ∑ ξ∈n exp [−β (H(ξ) − F)] ≡ exp [−β (Fn − F)] , (2.4) 23 2 Basics where we define Fn, the free energy of state n, in the last equality. The internal energy is defined as Un = ∑ ξ∈n P(ξ |n) = ∂ ∂β (βFn) , (2.5) where P(ξ |n) = exp [β (H(ξ) − F)] Peq(n) = exp [−β (H(ξ) − Fn)] (2.6) is the conditional probability to observe a microstate ξ given the mesostate n. The intrinsic entropy is Sn = ∑ ξ∈n P(ξ |n) ln P(ξ |n) = β2 ∂ ∂β Fn = β(Un − Fn), (2.7) where the last term is the Shannon entropy of the conditional probability P(ξ |n). In equilibrium, the likelihood of observing a forward and backward trajectory is the same by definition, i.e. the probability P(m, ∆t |n, 0) of observing the system in state n at time t = 0 and then later in state m at time t is the same as P(t, ∆t |m, 0). Therefore, it follows that P(m, t |n, 0)Peq(n) = P(n, t |m, 0)Peq(m). (2.8) For a small interval ∆t and using Eq. (2.4), we obtain the local detailed balance relation Peq(m) Peq(n) = wnm wmn = exp [−β (Fm − Fn)] = exp [−β (Em − En) + (Sm − Sn)] , (2.9) where we define the transition rate from state n to m as wnm ≡ lim ∆t→0 P(m, t |n, 0) ∆t . (2.10) We want to describe the dynamics of n(t) as a stochastic process. For that purpose, we will further assume that the dynamics of the microstates within each mesostate is much faster than the dynamics between the mesostates. In other words, we assume that the distribution of microstate within a mesostate relaxes into a local equilibrium distribution on a timescale much faster than the timescale for transitions between the mesostates. As a consequence, the dynamics of the mesostate n(t) becomes memoryless, i.e. the probability for a transition from n to m does not depend on the time spent in state n and states visited before n. Under this assumption of timescale separation, the dynamics of the mesostates can be described by a Markov process with the transition rates defined in Eq. (2.10). The time evolution of the probability P(n, t) to find the system in state n is governed by the master equation, which reads ∂ ∂t P(n, t) = ∑ m [P(m, t)wmn − P(n, t)wnm] . (2.11) 24 2.1 Stochastic thermodynamics for discrete state spaces Assuming the the system is ergodic, the probability P(n, t) will converge in the long time limit to the equilibrium distribution Peq(n), Eq. (2.4), which is time invariant. Indeed, plugging Eq. (2.4) into Eq. (2.11) and using Eq. (2.9) leads to 0 = ∑ m [Peq(m)wmn − Peq(n)wnm] . (2.12) In equilibrium, every terms in the sums vanish, i.e. there is no net probability flow across any link nm. 2.1.2 Open systems in contact with a heat bath and chemical reservoirs Without external control, a closed system in contact with a heat bath will always relax towards equilibrium. In contrast, open systems can achieve a nonequilibrium steady- state (NESS) by extracting energy from their environment. For example, biological systems rely on energy carriers such as adenosine triphosphates (ATP) and extract their energy by hydrolyzing these molecules. Here, we consider an open system coupled to chemical reservoirs. We group the main system and the chemical reservoirs into one large closed supersystem following the ideas from [58–60]. Chemical reservoirs are fully characterized by their chemical potential µα and we assume that they are much larger than the main system, i.e. we assume timescale separation between the dynamics of the reservoirs and the dynamics of the main system. From this assumption, the main system can be driven into a NESS by tuning the chemical potentials µα accordingly. The local detailed balance for the supersystem, Eq. (2.9), can be written as wnm wmn = exp [ −β ( Fm − Fn − ∑ α dαnmµα )] . (2.13) Here, we split the total free energy of the system into the free energy difference of the main system, (Fm − Fn), and the change in free energy of the reservoirs ∑ α dαnmµα. In a transition from n to m, molecules of species α at a chemical potential µα are supplied by the reservoirs (dnm < 0) or released in the reservoirs (dnm > 0). Analogously to Eq. (2.12), we can write a master equation for the probability distribu- tion P(n, t) of an open system as ∂ ∂t P(n, t) = ∑ m [P(m, t)wmn − P(n, t)wnm] = ∑ m jmn(t), (2.14) where we denote jnm(t) ≡ P(m, t)wmn − P(n, t)wnm (2.15) 25 2 Basics the net probability current from n to m. The Perron-Frobenius theorem [61] states that in the long time limit, P(n, t) will tend towards a unique stationary distribution Pss(n), which is obtained by solving the linear system of equations 0 = ∑ m [Pss(m)wmn − Pss(n)wnm] ≡ ∑ m jssmn , (2.16) where jssmn is the stationary probability current. In general, these currents do not vanish in a NESS. 2.1.3 Thermodynamics along an individual trajectory Stochastic thermodynamics is a framework that extends the concepts of classical ther- modynamics such as work, heat or entropy to the level of individual trajectories under nonequilibrium conditions [5–7]. In a closed system, Sekimoto [62] proposed that a change in internal energy along an individual trajectory must be compensated by a corresponding change in the energy of the heat bath. Specifically, the first law associated with a transition from n to m reads ∆Enm = −Qnm (2.17) where Qnm > 0 corresponds to heat dissipated in the heat bath. The total entropy change has three contributions. First, there is an increase in entropy of the bath ∆Sbath nm = βQnm (2.18) due to the heat Qnm dissipated. Second, there is a change of intrinsic entropy ∆Snm, defined in Eq. (2.7), due to a change of mesostate. The last contribution is the change in stochastic entropy [63] ∆Sstoch nm (t) ≡ ln P(n, t) P(m, t) . (2.19) If the probability distribution P(n, t) is not stationary, the stochastic entropy Sstoch(t) ≡ − ln P(n, t) will change even when the system remains in the same mesostate. These last two contributions can be combined as the change in system entropy ∆Ssys nm (t) ≡ ∆Snm + ∆Sstoch nm (t). (2.20) Consequently, the total entropy change for a transition from n to m can be written as as ∆Stot nm(t) = βQnm + ∆Ssys nm (t) = ln P(n, t)wnm P(m, t)wmn , (2.21) where the last equality follows from Eqs. (2.9) and (2.17). 26 2.2 Chemical reaction networks For an open system in contact with a heat bath and chemical reservoirs, the first law (2.17) along a transition from n to m becomes Enm = −Wchem mn −Qnm, (2.22) where we identify the chemical work performed on the reservoirs as Wchem mn ≡ ∑ α dαnmµα. (2.23) In a NESS, the total entropy production change along a stochastic trajectory n(t) can be written using Eq. (2.21) as ∆Stot [n(t)] = ∑ n,m Znm(t) ln P(n, t)wnm P(m, t)wmn , (2.24) where Znm(t) is a random variable that counts the number of transition from n to m in the time interval [0, t]. The mean rate of entropy production is [64] σ = ∑ mn P(n, t)wnm ln P(n, t)wnm P(m, t)wmn = ∑ n 0 and using the fact that (E−Y1 ) = − ( E+Y1 − 1 ) E−Y1 , Eq. (2.71) becomes∑ NY1 { ( E+Y1 − 1 ) [ k−1NY1 − k1E − Y1 ( NX1 − 2NY1 ) ( NX1 − 2NY1 − 1 )] + k2 ( E+ X1 E+ X1 E+X2 − 1 ) NX2 NY1 } δNY1,0 = 0. (2.74) The expressions proportional to k2 vanishes due to δNY1,0. Finally, we obtain the relation∑ NY1 NY1δNY1,0 = k1 k−1 NX1 ( NX1 − 1 ) , (2.75) Order 1 To solve the third line of Eq. (2.68), we note that the left hand side vanishes,∑ NY1 LfastP2(X, NY1, τ) = 0, (2.76) i.e. the fast operator Lfast only has a non-zero contribution for the lower order in ε in the perturbation serie, Eq. (2.67), the higher order terms evolve on the slow timescale τ. 37 2 Basics The right hand side is∑ NY1 [ ∂ ∂τ P0(X, NY1, τ) − L slowP1(X, NY1, τ) ] = 0. (2.77) The first term reduces to ∂/∂τP0(X, τ). Using Eq. (2.72), we can write the second term in the sum as∑ NY1 LslowP1(X, NY1, τ) = k1 ∑ NY1 ( E−Y1 − 1 ) ( NX1 − 2NY1 ) ( NX1 − 2NY1 − 1 ) P1(X, NY1, τ) = k1k2 ( E+ X1 E+ X1 E+X2 − 1 ) NX2 P0(X, τ) ∑ NY1 NY1δNY1,0. (2.78) Using Eq. (2.75), we obtain∑ NY1 LslowP1(X, NY1, τ) = k1k2 ( E+ X1 E+ X1 E+X2 − 1 ) NX2 NX1 ( NX1 − 1 ) P0(X, τ). (2.79) Combining these expressions, we obtain an effective CME ∂ ∂τ P0(X, τ) = k1k2 k−1 ( E+ X1 E+ X1 E+X2 − 1 ) NX2 NX1 ( NX1 − 1 ) P0(X, τ). (2.80) This equation is equivalent the chemical master Eq. (2.60) if we tune the rates such that k0 = k1k2/k−1 and in the limit ε → 0. We have thus shown that a trimolecular reaction of type (2.58) can be decomposed into two biomolecular reactions, Eq. (2.62), and have the same dynamics at the stochastic level assuming timescale separation. Note that our choice of biomolecular reactions is not unique. For a multimolecular reaction with arbitrary molecularity, Wilhelm proposed a systematic decomposition into a set of biomolecular reactions [82]. 2.2.5 Open chemical reaction networks Until now, we have dealt with closed chemical reaction networks where the total number of molecules ∑ i NXi is a conserved quantity. A closed system will always tend towards a unique equilibrium point in the long time limit. Open chemical reaction networks coupled to chemical reservoirs can achieve a nonequilibrium steady-state. We consider an open system with volume Ω consisting of a mixture of molecular intermediate species {X1, ..., XN } and chemostatted species {B1, ..., BC} supplied at con- stant concentrations by chemical reservoirs, also called chemostats. We consider that the reservoirs are ideal ones and provide a constant chemical potential µ j for all chemostatted 38 2.2 Chemical reaction networks species B j . We can connect the chemical potentials µ j with the concentrations b j using the equation for an ideal dilute solution, µ j = µ0 + β −1 ln b j b0 , (2.81) where µ0 is the chemical potential for standard conditions at concentration b0. A typical reaction can be written as N∑ i=1 si Xi + C∑ j=1 s j B j kρ −−−→ N∑ i=1 ri Xi + C∑ j=1 r j B j (2.82) The dynamical state of the system is given by the vector of intermediate species X(t) = ( NX1(t), ..., NXN (t) ) , where NXi (t) is the number of Xi molecules in the system at time t. Note that the state of the system does not depend explicitly on the concentrations of the chemostatted species. The time evolution of the system is described by the CME, Eq. (2.31). The propensity function associated with reaction (2.82) is given by aρ(X) = Ωkρ N∏ i NXi ! Ωsi ( NXi − si ) ! C∏ j=1 bsj j . (2.83) For thermodynamic consistency, the ratio of the reaction rate kρ and its inverse counterpart must fulfill the local detailed balance relation for open systems, Eq. (2.13). 39 3 First-order phase transitions in biochemical switches 3.1 Introduction In this chapter, we focus on biochemical switches that undergo a first-order phase transi- tion upon activation. From a deterministic perspective, this phase transition is associated with bistability where two stable steady states can coexist. In contrast, in the stochastic perspective, the steady-state is unique and is associated with a bimodal density distribu- tion [83, 84]. In the thermodynamic limit, the stochastic system will relax to the more stable fixed-point except at the bistable point [85, 86]. There are two timescales relevant for a biochemical switch: a fast relaxation to the nearest fixed-point and a slower transi- tion between the states, where the coarse-grained transition rates are proportional to the exponential of the inverse volume [87, 88]. In chemical reaction networks, bistability can occur only far from equilibrium since an equilibrium distribution will always be a Poisson distribution distribution, and, thus, have a single peak [77, 89]. The behavior of the entropy production at first- and second-order phase transitions has been investigated in many systems such as chemical networks or nonequilibrium Ising models [84–86, 90–100]. At first-order phase transitions, the entropy production rate has a discontinuity with respect to the thermodynamic force whereas at second-order phase transitions its first derivative has a discontinuity. Recently, it has been shown that the critical fluctuations of the entropy production diverge with a power-law with the volume at a second-order phase transition [98]. For first-order phase transitions, the behavior of entropy production fluctuations has not been investigated yet to the best of our knowledge. We will show that the fluctuations of the entropy production have an exponential volume-dependence at first-order phase transitions in chemical reaction networks. Our results are obtained for Schlögl’s model. First, we compute the entropy fluctuations numerically from the chemical master equation using standard large deviation techniques [101, 102]. Second, we compute the current fluctuations for a coarse-grained two-state model and show that the diffusion coefficient diverges at the bistable point with an exponential prefactor given by the height of the effective potential barrier separating the two fixed-points. As an alternative approach, we derive an analytical expression for the diffusion coefficient starting with the Fokker-Planck equation. The chapter is organized as follows. In Section 3.2, we introduce Schlögl’s model 41 3 First-order phase transitions in biochemical switches and define the entropy production. In Section 3.3, we consider the chemical master equation and compute the diffusion coefficient numerically. In Section 3.4, we introduce an effective two-state model and compute an analytical expression for the diffusion coefficient. In Section 3.5, we compute the diffusion coefficient starting with the Fokker- Planck equation. We conclude in Section 3.6. 3.2 Schlögl’s model and entropy production 3.2.1 Model definition The Schlögl model is a paradigmatic model for biochemical switches [51, 84]. It consists of a chemical species X in a volume Ω. The external bath contains two chemical species A and B at fixed concentrations a and b, respectively. The set of chemical reactions is 2X + A k1 −−⇀↽−− k−1 3X, B k2 −−⇀↽−− k−2 X, (3.1) where k1, k−1, k2 and k−2 are transition rates. The system is driven out of equilibrium due to a difference of chemical potential between A and B, which is written as ∆µ ≡ µA− µB. A cycle in which an X molecule is created with rate k1, and then degraded with rate k−2 leads to the consumption of a substrate A and generation of a product B. The thermodynamic force associated with this cycle is ∆µ ≡ ln k−2k1a k−1k2b , (3.2) where the temperatureT and Boltzmann’s constant kB are set to 1 throughout this chapter. 3.2.2 Entropy production Along a stochastic trajectory n(t), where n labels the state with n molecules of species X , the entropy production change of the heat bath can be identified as [6] ∆sbath = ZB(t) ln k−2 k2b + ZA(t) ln k1a k−1 . (3.3) Here, ZA(t) and ZB(t) are random variables which count transitions in the A and B channel, respectively. For example, ZB(t) increases by one if a B is produced, which 42 3.2 Schlögl’s model and entropy production happens if a reaction with rate k−2 takes place. Likewise, it decreases by one if a B is consumed, which happens if a reaction with rate k2 takes place, i.e. 2X + A k1 −−−→ 3X ZA → ZA + 1 , 3X k−1 −−−→ 2X + A ZA → ZA − 1 , X k−2 −−−→ B ZB → ZB + 1 , B k2 −−−→ X ZB → ZB − 1 . (3.4) In this chapter, we focus on ∆sbath, the extensive part of the total entropy production ∆stot ≡ ∆sbath + ∆s. The remaining part is the change in stochastic entropy [63] ∆s = − ln pnt (t) + ln pn0(0), (3.5) where pnt (t) is the probability to find the system in state nt at time t. Using Eqs. (3.3) and (3.4), we can write the mean entropy production rate per volume in the steady-state as [6] σ ≡ lim t→∞ 〈∆sbath〉 tΩ = lim t→∞ [ 〈ZB(t)〉 tΩ ln k−2 k2b + 〈ZA(t)〉 tΩ ln k1a k−1 ] = JB∆µ. (3.6) In the steady-state, the mean flux density of B molecules JB ≡ lim t→∞ 〈ZB(t)〉 tΩ = lim t→∞ 〈ZA(t)〉 tΩ (3.7) is equal to the flux of A molecules consumed. We can quantify the fluctuations of B molecules with the diffusion coefficient DB ≡ lim t→∞ 〈ZB(t)2〉 − 〈ZB(t)〉2 2tΩ = DA ≡ lim t→∞ 〈ZA(t)2〉 − 〈ZA(t)〉2 2tΩ , (3.8) where we prove the second equality, DB = DA, in Appendix 3.A. Specifically, we show there that ZA(t) and ZB(t) have the same cumulants. 43 3 First-order phase transitions in biochemical switches Finally, using Eqs. (3.3) and (3.8) we obtain the diffusion coefficient associated with the entropy production in the heat bath as Dσ ≡ lim t→∞ 〈(∆sbath)2〉 − 〈∆sbath〉2 2tΩ = DB ( ln k−2 k2b )2 + DA ( ln k1a k−1 )2 + lim t→∞ [ 〈ZA(t)ZB(t)〉 − 〈ZA(t)〉〈ZB(t)〉 2tΩ ] 2 ( ln k−2 k2b ) ( ln k1a k−1 ) = DB∆µ 2. (3.9) Since the stochastic entropy production is not extensive in time, this diffusion coefficient is equal to the one for the total entropy production. 3.3 Chemical master equation 3.3.1 Stationary solution The state of the system is fully determined by the total number n of X molecules. The time evolution of P(n, t), which is the probability to find the system in state n at time t, is governed by the chemical master equation (CME) ∂ ∂t P(n, t) = fn−1P(n − 1, t) + gn+1P(n + 1, t) − ( fn + gn)P(n, t) ≡ L0P(n, t). (3.10) Here, we define the rate parameters as fn = α+n + β + n ≡ ak1n(n − 1) Ω + bk2Ω, gn = α − n + β − n ≡ k−1n(n − 1)(n − 2) Ω2 + k−2n. (3.11) 44 3.3 Chemical master equation The system can reach a nonequilibrium steady-state with a distribution written as Pn. The analytical solution for this steady-state probability distribution reads Pn P0 = n−1∏ i=0 fi gi+1 , = f0 g1 √ f1g1 fngn exp [ 1 2 ln ( f1 g1 ) + n−1∑ i=2 ln ( fi gi ) + 1 2 ln ( fn gn )] (n ≥ 3), (3.12) where the normalization is given by P0 = 1 − ∞∑ j=1 Pj . (3.13) For a large number of states (Ω→∞), we can write x = n/Ω as a continuous variable and approximate the exponential by an integral using the trapezium rule, which is valid if fn/gn is bounded. Without loss of generality, we will choose parameters such that the fixed-points are far enough from the boundary, which means that the probability to find the system close to x = 0 is negligible. Note that the case where the solution does not vanish at the boundary is discussed in details in [88]. From Eq. (3.12), the continuous steady-state distribution can be written as p(x) ∝ exp [ −Ω ( φ0(x) + 1 Ω φ1(x) )] (3.14) where we define the non-equilibrium potential φ0(x) ≡ − ∫ x 0 dy ln ( f (y) g(y) ) , (3.15) and φ1(x) ≡ − 1 2 ln 1 f (x)g(x) . (3.16) Here, we have defined the total transition rates f (x) ≡ α+(x) + β+(x) = ak1x2 + k2b, g(x) ≡ α−(x) + β−(x) = k−1x3 + k−2x. (3.17) 45 3 First-order phase transitions in biochemical switches In the deterministic limit (Ω → ∞), we obtain the equation of the time evolution of the density, x̄ ≡ ∑ n nP(n, t)/Ω (3.18) as dx̄ dt = f (x̄) − g(x̄) (3.19) from the chemical master Eq. (3.10). In the steady-state, this equation has three solutions (x−, x0, x+). Bistability occurs when all solutions are real, we order the fixed-points as follows: 0 < x− < x0 < x+, where x± are stable (i.e. f ′(x±) < g′(x±)) and x0 is unstable (i.e. f ′(x0) > g′(x0)). 3.3.2 Behavior of the entropy production at the phase transition The mean entropy production rate, defined in Eq. (3.6), can be written using Eq. (3.17) as σ = ∆µ ∫ ∞ 0 dx [ β−(x) − β+(x) ] p(x), (3.20) In the thermodynamic limit, the stochastic systemwill relax to the more stable fixed-point except at the bistable point [85, 86]. Consequently, the rate of entropy production will be discontinuous with respect to the thermodynamic force at the bistable point. Specifically, with increasing ∆µ, σ will jump from [β−(x−) − β+(x−)] ∆µ to [β−(x+) − β+(x+)] ∆µ, which are the rates of entropy production at the two fixed-points x− and x+, respectively. We now derive an expression for the fluctuations of the entropy production. We follow an approach based on large deviation theory [101–103] and considered for Brownian ratchets in [104]. We want to compute the cumulants related to the number of produced B molecules. They are obtained through the scaled cumulant generating function (SCGF) α(λ) ≡ lim t→∞ 1 t ln〈eλZB〉 (3.21) where ZB is the time-integrated current of B molecules defined in Eq. (3.4). Note that α(λ) is unrelated to the transition rate defined in Eq. (3.11). From now on, we will drop the t dependence on ZB for readability. Expanding the generating function yields α(λ) = ΩJBλ +ΩDBλ 2 + O(λ3) (3.22) The SCGF can be obtained by considering the moment generating function g(λ, t) ≡ 〈eZB〉 = ∑ n g(λ, n, t)P(n, t) (3.23) 46 3.3 Chemical master equation where we define g(λ, n, t) ≡ 〈eλZB |n(t) = n〉 = ∑ ZB eλZBP(n, ZB, t), (3.24) which is conditioned on the final state of the trajectory n(t). The time evolution of this quantity is given by ∂ ∂t g(λ, n, t) = ∑ m Lnm(λ)g(λ,m, t). (3.25) where L(λ) is the tilted operator. Note that for λ = 0, it is identical to the operator generating the time evolution of the probability distribution in Eq. (3.10) We want to specify L(λ) for the chemical master Eq. (3.10). First, we write the time evolution of the probability distribution as P(n, ZB, t) = ∑ m wmnP(m, ZB − dmn, t) − wnmP(n, ZB, t) (3.26) wherewmn is the transition rate from statem to state n and dmn is the distancematrix which characterize how ZB changes during a transition. In the case of the CME, Eq. (3.26) reduces to P(n, ZB, t) = α−n+1P(n + 1, ZB, t) + α−n−1P(n − 1, ZB, t) + β−n+1P(n + 1, ZB + 1, t) + β+n−1P(n − 1, ZB − 1, t) − ( α+n + β + n + α − n + β − n ) P(n, ZB, t) (3.27) where the rates α±n and β±n are given by Eq. (3.11). Using Eq. (3.24), we can write ∂ ∂t g(λ, n, t) = ∑ m wmn ∑ ZB eλZBP(m, ZB − dmn, t) − wnmg(λ, n, t). (3.28) With a change of variable, we obtain ∂ ∂t g(λ, n, t) = ∑ m wmn ∑ YB eλ(YB+dmn)P(m,YB, t) − wnmg(λ, n, t) = ∑ m eλdmnwmng(λ,m, t) − wnmg(λ, n, t), (3.29) which specifies the tilted operator L(λ) defined in Eq. (3.25). 47 3 First-order phase transitions in biochemical switches The cumulants can be obtained by solving the eigenvalue equation L(λ)Q(λ) = α(λ)Q(λ), where L(λ) = L0 + L1λ + L2λ 2 + O(λ3) (3.30) and the distribution Q(λ) = P +Q1λ +Q2λ 2 + O(λ3). (3.31) Sorting by orders of λ, we obtain L0P = 0, L0Q1 + L1P = ΩJBP, L1Q1 + L0Q2 + L2P = ΩDBP +ΩJBQ1. (3.32) We multiply these equations with 〈1| on the left hand side and note that P is normalized, i.e. 〈1 | P〉 = 1, where 〈·|·〉 denotes the standard scalar product. We can compute the mean flux density JB = 1 Ω 〈1| L1 |P〉 = 1 Ω ∑ n ( β−n − β + n ) Pn (3.33) and the diffusion coefficient DB = 1 Ω ( 〈1| L1 |Q1〉 + 〈1| L2 |P〉 ) − JB 〈1 |Q1〉 = 1 Ω ∑ n ( ( β−n − β + n ) (Q1)n + 1 2 ( β−n + β + n ) Pn ) − JB ∑ n (Q1)n. (3.34) The mean rate of entropy production σ as well as its associated diffusion coefficient Dσ can be evaluated using Eqs. (3.6) and (3.9). 3.3.3 Numerical results Throughout this paper, we set the parameters to k1 = 1, k2 = 0.2, k−2 = 1, a = 1 and b = 1. The transition rate k−1 is computed from ∆µ and the generalized detailed balance relation Eq. (3.2), where ∆µ is a control parameter of the phase transition. In Fig. 3.1(a), we plot the stationary distribution p(x) which is bimodal in the vicinity of the phase transition (∆µbi ' 3.045). In Fig. 3.1(b), we plot the entropy production rate σ as a function of ∆µ. With increasing system size, σ gets steeper at the bistable point. In Fig. 3.1(c), we show that the first derivative of σ follows a power-law with an effective prefactor close to 1. In the thermodynamic limit (Ω → ∞), the entropy production rate becomes discontinuous as shown by Ge and Qian [85, 86]. 48 3.4 Two-state model The diffusion coefficient Dσ reaches a maximum at the bistable point and has an exponential volume-dependence. In Fig. 3.2(a), we compare simulations of the chemical master Eq. (3.10) using Gillespie’s algorithm [70] for three increasing sampling times T with the numerical solution obtained by solving the linear system given by Eq. (3.34). The systematic difference between these two methods is due to a limited sampling time. In the next section, we present an effective two-state model and derive an analytical expression for the diffusion coefficient. 3.4 Two-state model 3.4.1 Stationary solution In the bistable regime, the system has two timescales. First, it will relax towards the nearest stable fixed-points x± and fluctuate around it. Close to the fixed-points, the system can be modeled by stationary Gaussian processes. Specifically, the distribution p(x) can be expanded around its stable fixed-points x± as [84], p(x) ≈ ∑ x∗=(x−,x+) e−Ωφ0(x∗) ZG A2(x∗) exp [ −Ωφ′′0 (x∗) (x − x∗)2 2 ] , (3.35) where ZG ≡ ∑ x∗=(x−,x+) e−Ωφ0(x∗) √ 2π A2(x∗) √��φ′′0 (x∗)��Ω . (3.36) Second, as stochastic fluctuations are always present, the system will at some point in time reach the unstable fixed-point x0 beyond which it can relax towards the other fixed- point. Based on this behavior, the infinite-state system Eq. (3.10) can be coarse-grained into a two-state process between the stable fixed-points x± [105]. The transition rates from x± to x∓ depend exponentially on the system size and are given explicitly by [87, 88] r± = e−Ω[φ0(x0)−φ0(x±)] f (x±) √ −φ′′0 (x0)φ ′′ 0 (x±) 2πΩ . (3.37) In Appendix 3.B, we derive these transition rates following Hinch and Chapman [88]. 3.4.2 Behavior of the entropy production at the phase transition Wewant to characterize the fluctuations of the entropy production for the two-statemodel. We consider the thermodynamic flux JB and its associated diffusion coefficient DB defined 49 3 First-order phase transitions in biochemical switches 0 1 2 3 4 5 0 2 4 6 x p(x) ∆µ = 3.00 ∆µ = 3.05 ∆µ = 3.10 ∆µ = 3.15 (a) 3 3.02 3.04 3.06 3.08 3.1 3.12 3.14 0 2 4 6 8 10 ∆µ σ Ω = 50 Ω = 100 Ω = 200 Ω = 1000 (b) 200 400 600 800 1,000 500 1,000 1,500 Ω ∂σ ∂∆µ ∝ Ω0.97±0.01 (c) Figure 3.1: Phase transition in the Schlögl model. (a) Stationary distribution of chemical species X for Ω = 100 and different values of ∆µ. (b) Mean entropy produc- tion rate σ as a function of ∆µ for different system sizes Ω. (c) Maximum of the first derivative of σ as a function of the system size Ω. Parameters are given in the main text. 50 3.4 Two-state model 3.05 3.055 3.06 3.065 3.07 3.075 3.08 15 20 25 ∆µ lnDσ Gillespie (T = 104) Gillespie (T = 105) Gillespie (T = 106) CME 3.07 3.08 18.6 18.8 19 (a) 2.95 3 3.05 3.1 3.15 3.2 0 5 10 15 20 ∆µ lnDσ CME Two-state FP (b) 200 400 600 800 1,000 20 40 60 80 100 Ω lnDbi σ CME Two-state FP (c) Figure 3.2: Behavior of the diffusion coefficient Dσ close to the bistable point. (a) Diffusion coefficient from simulations using Gillespie’s algorithm for 104 trajectories with a sampling time ofT = 104, 105, 106. For the CME, we solve the linear system given by Eq. (3.34) numerically for Ω = 100. (b) Diffusion coefficient for the two-state model obtained by evaluating Eq. (3.43) and for the Fokker-Planck equation (FP) by evaluating Eq. (3.74) (c) Finite-size scaling of the maximum of the diffusion coefficient Dbi σ . The corresponding slopes are given in Table 3.1. 51 3 First-order phase transitions in biochemical switches in Eqs. (3.7) and (3.8). There are two contributions to both of these quantities. First, the system can fluctuate around one of the fixed-point, which is modeled by a Gaussian process Eq. (3.35). Within the state x±, the thermodynamic flux J± and its diffusion coefficient D± can be computed exactly [102]. Second, the system can jump from state x± to x∓ on the largest timescale. For simplicity, we will assume that a transition from x± to x∓ produces an average flux of B± molecules, where we expect B± ∼ O(Ω). The probability to be in state x± is p± ≡ r∓ r− + r+ . (3.38) We will compute DB using large deviation theory as introduced in Section 3.3.2. In Appendix 3.C, we compute DB without relying on large deviation theory. Combining these two contributions, the tilted operator for this two-system system reads [106] L(λ) = ( −r− + J−λ + D−λ2 r+eλB+ r−eλB− −r+ + J+λ + D+λ2 ) . (3.39) The maximal eigenvalue of L(λ) is α(λ) = TrL(λ)/2 + √ (TrL(λ))2 /4 − DetL(λ), (3.40) where Tr and Det denote the trace and the determinant, respectively. From Eq. (3.22), the average flux of B is given by JB = ∂α(λ) ∂λ ����� λ=0 = p−J− + p+J+ + r−p−(B− + B+) = p−J− + p+J+ + O ( e−Ω|∆φ| ) (3.41) where ∆φ0 = φ0(x0) − φ0(x−). The diffusion coefficient reads DB = 1 2 ∂2α(λ) ∂λ2 ����� λ=0 = p−p+ (J− − J+)2 r− + r+ + p−D− + p+D+ + p−p+ (B− + B+) (p+ − p−)(J− − J+) + 1 2 (B− + B+) 2 p−p+ (r−p+ + p−r+) = p−p+ (J− − J+)2 r− + r+ + O(Ω). (3.42) 52 3.4 Two-state model Method Exponential prefactor δ (Dσ ∝ eδΩ) Eq. (3.34) 0.0846 ± 0.0005 Eq. (3.43) 0.0846 ± 0.0005 Eq. (3.74) 0.0845 ± 0.0006 Eq. (3.44) 0.0823 Table 3.1: Scaling of the diffusion coefficient Dσ obtained with the CME, Eq. (3.34), the two-state model, Eq. (3.43), the Fokker-Planck equation, Eq. (3.74) and by evaluating the height of the effective potential barrier separating the two fixed-points, Eq. (3.44). The maximum logarithm of the diffusion coefficient is fitted with an linear function of the system size Ω and the errors is given with 95% confidence bounds. We insert the transition rates Eq. (3.37) into the previous expression and obtain DB �� ∆µbi = p−p+ (J− − J+)2 eΩ[φ0(x0)−φ0(x±)]πΩ f (x−) √ −φ′′0 (x0)φ ′′ 0 (x±) , (3.43) where f (x−) = f (x+) and φ0(x−) = φ0(x+) at the bistable point. As a main result, we have thus shown that the diffusion coefficient scales as DB �� ∆µbi ∝ eΩ[φ0(x0)−φ0(x±)], (3.44) where the exponential prefactor is the height of the effective potential barrier between the two fixed-points. The mean rate of entropy production σ as well as its associated diffusion coefficient Dσ can be evaluated using Eqs. (3.6) and (3.9). Here, we have assumed that the average flux of B molecules B± produced during a jump from x± to x∓ is known and inserted these values in Eq. (3.39). In fact, a trajectory from x− to x+ will produce a path-dependent flux of B±. When we perform a coarse- graining of the CME into a two-state model, we lose this information. Nevertheless, as we are only interested in the leading terms of DB, we have shown that the contributions from jumps between the fixed-points B± can be neglected close to the bistable point for large system sizes. 3.4.3 Numerical results We compare the analytical results with numerical evaluations the CME. In Fig. 3.2(b), we show Dσ, Eq. (3.43), and compare it with the numerical results from the CME. We find that Dσ evaluated for the two-state model almost matches the CME close to the bistable point. In Fig. 3.2(c), we show that the diffusion coefficient has an exponential volume- dependence at the bistable point. In Table 3.1, we compare the scaling of the maximum 53 3 First-order phase transitions in biochemical switches of the diffusion coefficient obtained numerically and by evaluating Eq. (3.44). The difference between the numerical prefactors and our analytical expression, Eqs. (3.43) and (3.44), is due to finite-size effects. 3.5 Fokker-Planck equation 3.5.1 Stationary solution and diffusion dilemma The Fokker-Planck equation is usually obtained by truncating the Kramers-Moyal ex- pansion of the chemical master Eq. (3.10). Retaining only the first two terms leads to ∂ ∂t p(x, t) = − ∂ ∂x A1(x)p(x, t) + 1 2Ω ∂2 ∂x2 A2(x)p(x, t) (3.45) where A1(x) = f (x) − g(x) and A2(x) = f (x) + g(x) (3.46) are the first two Kramers-Moyal moments. The stationary distribution of this Fokker- Planck equation is p(x) ∝ exp [ −ΩΨ0(x) ] , (3.47) where Ψ0(x) = −2 ∫ x 0 dy A1(y) A2(y) . (3.48) However, this potential Ψ0(x) differs from the CME potential φ0(x) in leading Ω orders. As discussed in Section 2.2.3, it is not a suitable choice to describe the transitions between the fixed-points, which are rare events for large system sizes. Hänggi et al. [87] proposed a different diffusion process described by a Fokker-Planck equation of the form ∂ ∂t p(x, t) = − ∂ ∂x [ L(x)χ0(x)p(x, t) ] + 1 Ω ∂2 ∂x2 [ L(x)p(x, t) ] ≡ LFP 0 p(x, t), (3.49) where we define, using Eq. (3.15), the generalized thermodynamic force χ0(x) ≡ − ∂φ0(x) ∂x (3.50) and an effective diffusion coefficient expressed with a transport coefficient [107, 108] L(x) ≡ f (x) − g(x) ln f (x) − ln g(x) . (3.51) The stationary solution of Eq. (3.49) is p(x) ∝ exp [−Ωφ0(x)] , (3.52) 54 3.5 Fokker-Planck equation which is equivalent to Eq. (3.14) in leading Ω orders. φ0(x) is the large deviation rate function of p(x) [86, 102, 103]. Thus, in the thermodynamic limit (Ω → ∞), Eq. (3.49) has the same large deviation function as the CME. In other words, Eq. (3.49) will reproduce the correct long-term asymptotic dynamics. However, Eq. (3.49) will not correctly reproduce the stochastic dynamics on the shorter timescale. As shown in Section 2.2.3, the correct diffusion coefficient on the shorter timescale is given by the second Kramers-Moyal moment, which differs from Eq. (3.51) at leading order in Ω. Therefore, there is dilemma when choosing the diffusion term in the Fokker-Planck equation [109, 110]. It is not possible to approximate the CME with a diffusion process that describes correctly the short-time stochastic dynamics and the long-time global dynamics. As shown in Section 3.4, the leading terms of the diffusion coefficient DB are due to transitions between the fixed-points, thus, Eq. (3.49) is the appropriate Fokker- Planck equation to characterize DB close to the bistable point. 3.5.2 Behavior of the entropy production at the phase transition In Section 3.3.2, we derived the tilted operator L(λ) for the CME, Eq. (3.29). Analo- gously, one can construct a tilted operator L(λ) for the Fokker-Planck equation, see [103, 111] for additional details. For the time-integrated current ZB, the tiled operator is given by LKM(λ) = − ( ∂ ∂x + λΩζB ) A1(x) + 1 2Ω ( ∂ ∂x + λΩζB )2 A2(x) (3.53) where ζB is an operator that only selects transitions in the B channel occurs. For simplicity, we first consider the Kramer-Moyal Fokker-Planck equation. ζB acts on the Fokker-Planck coefficients as follows: ζB A1(x) ≡ ( β−(x) − β+(x) ) , ζB A2(x) ≡ ( β+(x) + β−(x) ) . (3.54) The tilted operator coefficient LKM 0 is given by Eq. (3.45), LKM 1 ≡ Ω [ β−(x) − β+(x) ] + ∂ ∂x ( β+(x) + β−(x) ) (3.55) and LKM 2 = Ω 2 [ β+(x) + β−(x) ] . (3.56) In Appendix 3.D, we derive Eqs. (3.55) and (3.56) by discretizing the tilted operator of the CME. 55 3 First-order phase transitions in biochemical switches We now consider the Fokker-Planck Eq. (3.49), where the tilted operator reads LFP(λ) = − ( ∂ ∂x + λΩζB ) L(x)χ0(x) + 1 2Ω ( ∂ ∂x + λΩζB )2 2L(x) (3.57) In order to write the first-order coefficient LFP 1 , we must first compute how the operator ζB acts on the diffusion coefficient, 2L(x). The drift term reads A1(x) = L(x)χ0(x). (3.58) This expression can also be expressd as [112] A1(x) = 1 2P(x) ∞∑ n=1 (−1)n+1 n!Ωn ∂n ∂xn An+1(x)P(x) (3.59) where the Kramers-Moyal moments are given by An(x) = f (x) + (−1)n g(x) (3.60) and P(x) = Ωp(Ωx) (3.61) By expanding the right hand side of Eq. (3.59) in powers of Ω−1, we can express the transport coefficient L(x) as a sum of the Kramers-Moyal moments using Eq. (3.58) leading to L(x) = 1 2 ∞∑ n=0 1 (n + 1)! An+2(x) (χ0(x))n . (3.62) Using Eq. (3.54), we can compute how the operator ζB acts on L(x). We define the transport coefficient for the B channel as LB(x) ≡ ζBL(x), = 1 2 ∞∑ n=0 1 (n + 1)! ( β+(x) + (−1)n β−(x) ) (χ0(x))n , = L(x) 2 (β−(x) f (x) + β+(x)g(x)) f (x)g(x) . (3.63) The tilted operator coefficient LFP 0 is given by Eq. (3.45), LFP 1 ≡ Ω [ β−(x) − β+(x) ] + 2 ∂ ∂x LB(x) (3.64) 56 3.5 Fokker-Planck equation and LFP 2 ≡ ΩLB(x). (3.65) We can write the diffusion coefficient, defined for the CME in Eq. (3.34), as DFP B = ∫ ∞ 0 dx { [ β−(x) − β+(x) ] q1(x) + 1 2 [ β−(x) − β+(x) ] p(x) } (3.66) where the function q1(x) is obtained by solving the differential equation LFP 0 q1(x) = ( JB − L FP 1 ) p(x). (3.67) Writing the solution as q1(x) = C(x)p(x), it can be shown (e.g. using Mathematica) that the solution is ∂C(x) ∂x = C1 + ∫ x 0 dy p(y)L(y)F(y) p(x)L(x) (3.68) where C1 = 0 due to the boundary condition (x = 0) and F(y) = { Ω [ JB f (x)g(x) + β+(x)g(x)2 − β−(x) f (x)2 f (y)g(y)L(y) ] − ∂ ∂y [ β+(y) f (y) + β−(y) g(y) ] } . (3.69) Close to the bistable point, the system has three fixed-points, which we order as follows: 0 < x− < x0 < x+, where x± are stable (i.e. f ′(x±) < g′(x±)) and x0 is unstable (i.e. f ′(x0) > g′(x0)). In the limit of large system sizes, we can write p(x) as a sum of two stationary Gaussian processes, Eq. (3.35). Using a saddle-point approximation, Eq. (3.68) becomes ∂C(x) ∂x ≈  ∑ x∗=(x−,x+) Θ (x − x∗) F(x∗)e−Ωφ0(x∗) √ 2π��φ′′0 (x∗)��Ω  eΩφ0(x0) √ 2π��φ′′0 (x0) ��Ω exp [ −Ωφ′′0 (x0) (x − x0) 2 2 ] . (3.70) where Θ(x) is the Heaviside function. We can integrate this expression and obtain C(x) ≈ Θ (x − x−) eΩ[φ0(x0)−φ0(x−)] 2πF(x−) Ω √��φ′′0 (x−)φ′′0 (x0) �� + C0, (3.71) where C0 is a normalization constant. By multiplying Eq. (3.71) with Eq. (3.35), we obtain q1(x) ≈ 1 ZG 2πF(x−)eΩ[φ0(x0)−φ0(x−)−φ0(x+)] Ω f (x+) √��φ′′0 (x−)φ′′0 (x0) �� exp [ −Ωφ′′0 (x+) (x − x+)2 2 ] +C0p(x). (3.72) 57 3 First-order phase transitions in biochemical switches where the normalization C0 is chosen such that∫ ∞ 0 dx q1(x) = 0. (3.73) Close to the bistable point, we show that the diffusion coefficient diverges as an exponen- tial function of the system size. We can neglect the term ∫ ∞ 0 dx L2p(x) as it is of order 1 as shown in Section 3.4.2. Eq. (3.66) then becomes DFP B ≈ ∫ ∞ 0 dx [ β−(x) − β+(x) − JB ] q1(x) ≈ 1 ZG ( 2π Ω )3/2 F(x−)eΩ[φ0(x0)−φ0(x−)−φ0(x+)] f (x+) √��φ′′0 (x−)φ′′0 (x0)φ ′′ 0 (x+) �� [ β−(x+) − β+(x+) − JB ] . (3.74) When the system is bistable, the diffusion coefficient diverges as an exponential law of the system size, namely, DFP B ∝ eΩ[φ0(x0)−φ0(x−)−φ0(x+)−Ω−1 lnZG], (3.75) where the normalization constant ZG is given by Eq. (3.36). At the bistable point ∆µbi, both stable fixed-points have the same potential height [113], i.e. φ0(x−) = φ0(x+). The fourth term in Eq. (3.75) reduces to lim Ω→∞ Ω −1 lnZG ���� ∆µbi = −φ0(x−). (3.76) By inserting this expression into Eq. (3.75), we show that the diffusion coefficient scales as DFP B �� ∆µbi ∝ eΩ[φ0(x0)−φ0(x−)], (3.77) where the exponential prefactor is the same as in Eq. (3.44). 3.5.3 Numerical results In Fig. 3.2, we plot the diffusion coefficient obtained by evaluating Eq. (3.74) with the parameters given in Section 3.3.3. Close to the bistable point, Eq. (3.74) matches the results obtained for the two-state model, Eq. (3.43) Below the bistable point, Eq. (3.74) deviates from the CME and two-state model expressions. This could be explained by the use of saddle-point approximations, which are only valid when the distribution p(x) has two “distinct” peaks. 58 3.6 Conclusion 3.6 Conclusion We have investigated the fluctuations of the entropy production at the phase transition occurring in a paradigmatic model of biochemical switches. A control parameter for this phase transition is the thermodynamic force driving the system out of equilibrium. The mean entropy production rate has a discontinuity with respect to the therm