Thermodynamics of Chemical Waves

TL;DR

This paper develops a thermodynamic framework for reaction-diffusion systems to analyze chemical waves, defining a nonequilibrium free energy density and examining wave dynamics in biological contexts.

Contribution

It introduces a rigorous thermodynamic theory for chemical waves, including the definition of a nonequilibrium free energy and its balance equation, advancing understanding of biological information transfer.

Findings

Defined thermodynamic potential as nonequilibrium free energy density.

Analyzed wavefronts and oscillating waves using Fisher-Kolmogorov and Brusselator models.

Linked thermodynamics of chemical waves to biological information transfer.

Abstract

Chemical waves constitute a known class of dissipative structures emerging in reaction-diffusion systems. They play a crucial role in biology, spreading information rapidly to synchronize and coordinate biological events. We develop a rigorous thermodynamic theory of reaction-diffusion systems to characterize chemical waves. Our main result is the definition of the proper thermodynamic potential of the local dynamics as a nonequilibrium free energy density and establishing its balance equation. This enables us to identify the dynamics of the free energy, of the dissipation, and of the work spent to sustain the wave propagation. Two prototypical classes of chemical waves are examined. From a thermodynamic perspective, the first is sustained by relaxation towards equilibrium and the second by nonconservative forces generated by chemostats. We analytically study step-like waves, called wavefronts, using the Fisher-Kolmogorov equation as representative of the first class and oscillating waves in the Brusselator model as representative of the second. Given the fundamental role of chemical waves as message carriers in biosystems, our thermodynamic theory constitutes an important step toward an understanding of information transfers and processing in biology.

Tables

Table 1. Table 1: Set of units measures adopted for the Brusselator model

Physical quantity	Unit measures
time	$1/k_{+1}$
space	$\sqrt{D_1/k_{+1}}$
concentration	$k_{+1}/k_{+2}$

Equations

$${\partial }_t{Z}^{\alpha }={\mathbb{S}}_{\rho }^{\alpha }{j}^{\rho }-\nabla \cdot {\mathbf{J}}^{\alpha }+{\delta }_{{\alpha }_y}^{\alpha }{I}^{{\alpha }_y},$$

$$j^{\pm \rho }=k_{\pm \rho }\prod _{\alpha }(Z^{\alpha }{)}^{{\mathbb{S}}_{\pm \rho }^{\alpha }}.$$

$${\mathbf{J}}^{\alpha }=-D^{\alpha }\nabla Z^{\alpha },$$

$$l_{\alpha }^{\lambda }{\mathbb{S}}_{\rho }^{\alpha }=0\forall \rho \in \mathcal{R},$$

$$L^{{\lambda }_u}={\int }_V\mathrm{d}\mathbf{r}{l}_{{\alpha }_x}^{{\lambda }_u}{Z}^{{\alpha }_x}$$

$$M^{{\alpha }_{y_b}}:={\hat{l}}_{{\lambda }_b}^{{\alpha }_{y_b}}{\int }_V\mathrm{d}\mathbf{r}{l}_{\alpha }^{{\lambda }_b}{Z}^{\alpha }$$

$$k_{+\rho }/k_{-\rho }=\prod _{\alpha }(Z_{\text{eq}}^{\alpha }{)}^{{\mathbb{S}}_{\rho }^{\alpha }}.$$

$$\frac{k_{+\rho }}{k_{-\rho }}=\exp \left(-\frac{{\mathbb{S}}_{\rho }^{\alpha }{\mu }_{\alpha }^{\circ }}{RT}\right).$$

$${\mathrm{d}}_t\mathcal{G}=-T\dot{\Sigma }+{\dot{W}}_{\text{driv}}^{\mathcal{G}}+{\dot{W}}_{\text{nc}},$$

$$\mathcal{G}:=G-{\mu }_{{\alpha }_{y_b}}^{\text{ref}}{M}^{{\alpha }_{y_b}}$$

$$g=Z^{\alpha }{\mu }_{\alpha }-RT{Z}^{\mathcal{S}}$$

$${\dot{\Sigma }}_{\text{rct}}={\int }_V\mathrm{d}\mathbf{r}{\dot{\sigma }}_{\text{rct}}\ge 0,$$

$$T{\dot{\sigma }}_{\text{rct}}=-({\mu }_{\alpha }{\mathbb{S}}_{\rho }^{\alpha })j^{\rho }\ge 0,$$

$${\dot{W}}_{\text{driv}}^{\mathcal{G}}=-{\mathrm{d}}_t{\mu }_{{\alpha }_{y_b}}^{\text{ref}}{M}^{{\alpha }_{y_b}}.$$

$${\dot{W}}_{\text{nc}}={\int }_V\mathrm{d}\mathbf{r}{\dot{w}}_{\text{nc}},$$

$${\dot{w}}_{\text{nc}}={\mathcal{F}}_{{\alpha }_y}{I}^{{\alpha }_y},$$

$${\mathrm{d}}_t\mathcal{G}=-T\dot{\Sigma }\le 0.$$

$$\mathcal{G}={\mathcal{G}}_{\text{eq}}+RT{\int }_V\mathrm{d}\mathbf{r}\mathcal{L}(\{Z^{\alpha }(\mathbf{r},t)\}\parallel \{Z_{\text{eq}}^{\alpha }\})$$

$$\mathcal{L}(\{Z^{\alpha }\}\parallel \{Z_{\text{eq}}^{\alpha }\})=\sum _{\alpha \in \mathcal{S}}{Z}^{\alpha }\ln \left(\frac{Z^{\alpha }}{Z_{\text{eq}}^{\alpha }}\right)-\left(Z^{\alpha }-Z_{\text{eq}}^{\alpha }\right),$$

$$Z_{\alpha }^{\text{eq}}=\exp (({\mu }_{\alpha }^{\text{ref}}-{\mu }_{\alpha }^{\circ })/RT)$$

$$\text{mathpzc}g(\mathbf{r},t):=g(\mathbf{r},t)-{\mu }_{{\alpha }_{y_b}}^{\text{ref}}{\hat{l}}_{{\lambda }_b}^{{\alpha }_{y_b}}{l}_{\alpha }^{{\lambda }_b}{Z}^{\alpha }(\mathbf{r},t),$$

$${\partial }_t\text{mathpzc}g=-T\dot{\sigma }+{\dot{w}}_{\text{driv}}^{\text{mathpzc}g}+{\dot{w}}_{\text{nc}}-\nabla \cdot {\mathbf{J}}^{\text{mathpzc}g},$$

$${\dot{w}}_{\text{driv}}^{\text{mathpzc}g}=-{\mathrm{d}}_t{\mu }_{{\alpha }_{y_b}}^{\text{ref}}{\hat{l}}_{{\lambda }_b}^{{\alpha }_{y_b}}{l}_{\alpha }^{{\lambda }_b}{Z}^{\alpha }.$$

$${\mathbf{J}}^{\text{mathpzc}g}=({\mu }_{\alpha }-{\mu }_{{\alpha }_{y_b}}^{\text{ref}}{\hat{l}}_{{\lambda }_b}^{{\alpha }_{y_b}}{l}_{\alpha }^{{\lambda }_b}){\mathbf{J}}^{\alpha }$$

$$\text{mathbbm}g(\mathbf{r},t):=g(\mathbf{r},t)-{\mu }_{\alpha }^{\text{ref}}{Z}^{\alpha }(\mathbf{r},t),$$

$$\text{mathbbm}g(\mathbf{r},t)=\text{mathbbm}{g}_{\text{eq}}+RT\mathcal{L}(\{Z^{\alpha }(\mathbf{r},t)\}\parallel \{Z_{\text{eq}}^{\alpha }\})$$

$${\partial }_t\text{mathbbm}g=-T\dot{\sigma }+{\dot{w}}_{\text{driv}}^{\text{mathbbm}g}+{\dot{w}}_{\text{nc}}-\nabla \cdot {\mathbf{J}}^{\text{mathbbm}g}$$

$${\dot{w}}_{\text{driv}}^{\text{mathbbm}g}=-{\mathrm{d}}_t{\mu }_{\alpha }^{\text{ref}}{Z}^{\alpha }$$

$${\mathbf{J}}^{\text{mathbbm}g}=({\mu }_{\alpha }-{\mu }_{\alpha }^{\text{ref}}){\mathbf{J}}^{\alpha }$$

$$\text{mathbbm}G(t):={\int }_V\mathrm{d}\mathbf{r}\text{mathbbm}g(\mathbf{r},t),$$

Full text

I Introduction

Chemical waves, or traveling waves, are complex concentration patterns of chemical species moving in space with a constant velocity and without changes of shape Murray (2002). They are essential for communication in biosystems both at the intracellular and intercellular level, and play a crucial role for the synchronization and coordination of biological events. By using chemical reactions, chemical waves are able to spread signals more rapidly compared to simple diffusion Deneke and Di Talia (2018). The information they carry is encoded not only in the identity of the chemical species but also in other features of the concentration patterns such as the amplitude for step-like waves, called wavefronts, or the wavenumber for periodic waves Purvis and Lahav (2013); Behar and Hoffmann (2010). Calcium waves, for instance, are highly versatile intracellular messenger creating different spatiotemporal patterns Berridge et al. (2003, 2000) that regulates several cellular activities over a wide range of time scales Carafoli et al. (2001).

From a thermodynamic standpoint, chemical waves are out of equilibrium processes requiring a continual influx of energy and chemicals. They belong to the broader class of dissipative structures occurring in reaction-diffusion systems Prigogine and Nicolis (1971); Nicolis and Prigogine (1977). A lot of work has been dedicated to investigating their origin and their relative stability, but mainly from a dynamical point of view Cross and Hohenberg (1993); Ross (2008). However, a proper understanding of the energetic cost needed to support chemical waves is missing. We fill this gap in this work by developing a thermodynamic theory of chemical waves.

We start in Sec. II by providing a local formulation of thermodynamics of reaction-diffusion systems. Our theory is based on recent advances in the thermodynamic description of chemical reaction networks Rao and Esposito (2016, 2018a) and reaction-diffusion systems Falasco et al. (2018) inspired by stochastic thermodynamics Jarzynski (2011); den Broeck and Esposito (2015) and its links to information theory Parrondo et al. (2015). In Ref. Falasco et al. (2018), global thermodynamic quantities were defined in terms of the reaction-diffusion dynamics to study the energetics of dissipative structures. However, to charaterize localized processes such as chemical waves, a local formulation is needed in the spirit of the phenomenological nonequilibrium thermodynamics developed in the last century de Groot and Mazur (1984); Glansdorff and Prigogine (1971). The theory we develop does so but without relying on any linear approximation of the thermodynamic currents in the forces. It is systematically built on top of the dynamics and thus valid arbitrarily far from equilibrium. The thermodynamic quantities are local and expressed in terms of densities. A crucial step is the identification of the proper thermodynamic potential for the local dynamics, i.e., a free energy density. In a way which is reminiscent of equilibrium thermodynamics when passing from canonical to grand canonical ensembles, this potential is constructed from the Gibbs free energy density by eliminating the energetic contribution due to the local diffusion of the chemical species and by making use of the conservation laws of the chemical reaction network Polettini and Esposito (2014); Rao and Esposito (2018b, a). Its balance equation provides a local formulation of the second law and decomposes the evolution of the free energy density into a flow of free energy caused by diffusion and three source terms. The first two source terms are chemical works produced respectively by time dependent driving and by nonconservative forces while the third one is a sink term due to entropy production. For undriven detailed balanced systems (i.e., in absence of work), this free energy is minimized as the system relaxes to equilibrium.

We proceed in Sec. III by specializing our theory to study the thermodynamic cost of propagating and sustaining chemical waves. We first identify the conditions needed for an open reaction diffusion system to allow the existence of chemical waves. We then show in general that while the free energy is exclusively changing due to propagation, work and entropy production also quantify the cost for sustaining the waves. We subsequently consider in Sec. IV two prototypical wave patterns which we characterize analytically, i.e., wavefronts in the Fisher-Kolmogorov equation Fisher (1937); Murray (2002) and periodic/oscillating waves in the Brusselator model Prigogine and Lefever (1968); Auchmuty and Nicolis (1976). We find that when wavefronts are caused by relaxation towards equilibrium in detailed balanced systems, as Fisher-Kolmogorov waves, they propagate efficiently from an energetic standpoint because dissipation is solely localized at the front. Instead, oscillating waves, as in the Brusselator model, are significantly more energy costly as they need to be sustained by nonconservative forces (non detailed balanced systems) which cause dissipation throughout the system.

The implications of our work for information processing in biosystems as well as possible future developments are discussed in Sec. V.

II Reaction-Diffusion Systems

A Reaction-Diffusion System (RDS) is considered here as a dilute ideal mixture of chemical species $\ch{Z}_{\alpha}$ ($\alpha\in\mathcal{S}$) diffusing within a volume $V$ with impermeable boundaries $\partial V$ and undergoing elementary reactions $\rho\in\mathcal{R}$. Each chemical species $\alpha$ is classified either as a controlled/chemostatted species ($\alpha_{y}\in Y$) if it is exchanged with particular reservoirs called chemostats controlling its concentration or as internal species $(\alpha_{x}\in X)$ otherwise Rao and Esposito (2016). It is characterized by its concentration field $Z^{\alpha}(\boldsymbol{r},t)$ which follows the reaction-diffusion equation:

[TABLE]

where $\delta^{\alpha}_{\alpha_{y}}$ is the Kronecker delta and repeated upper-lower indices imply the summation over all the allowed values of the indices in all the paper. $j^{\rho}$ gives the concentration variation due to the chemical reaction $\rho$, $\mathbb{S}^{\alpha}_{+\rho}\ch{Z}_{\alpha}\mathbb{\rightleftharpoons}\mathbb{S}^{\alpha}_{-\rho}\ch{Z}_{\alpha}$. The stoichiometric coefficient $\mathbb{S}^{\alpha}_{+\rho}$ ($\mathbb{S}^{\alpha}_{-\rho}$) specifies the number of molecules $\alpha$ involved in the forward reaction $+\rho$ (backward reaction $-\rho$) whose chemical reaction current $j^{+\rho}$ ($j^{-\rho}$) satisfy the mass-action kinetics de Groot and Mazur (1984); Pekař (2005),

[TABLE]

The so-called stoichiometric matrix $\mathbb{S}^{\alpha}_{\rho}=\mathbb{S}^{\alpha}_{-\rho}-\mathbb{S}^{\alpha}_{+\rho}$ gives the full variation of the number of molecules $\alpha$ upon the reaction $\rho$, while $j^{\rho}=j^{+\rho}-j^{-\rho}$ specifies the net reaction current. The Fick’s diffusion currents

[TABLE]

with the diffusion coefficient $D^{\alpha}$ for the species $\alpha$, describe the transport of chemical species within the volume $V$. They vanish at the boundaries $\partial V$, i.e., $\int_{V}\mathrm{d}r\text{ }\nabla\cdot\boldsymbol{J}^{\alpha}=0$ $\forall\alpha$. If the RDS is open, the controlled species $\alpha_{y}\in Y$ are exchanged with the chemostats according to the external currents $I^{\alpha_{y}}(\boldsymbol{r},t)$. These latter vanish instead in a closed RDS.

The left null eigenvectors of the stoichiometric matrix,

[TABLE]

are the so-called conservation laws $l^{\lambda}$ Polettini and Esposito (2014); Rao and Esposito (2016) defining the quantities $L^{\lambda}=\int_{V}\mathrm{d}\boldsymbol{r}\text{ }l^{\lambda}_{\alpha}Z^{\alpha}$ that are conserved if the RDS is closed Falasco et al. (2018): $\mathrm{d}_{t}L^{\lambda}=0$. When the system is open, the set of conservation laws $l^{\lambda}$ is split into two disjoint subsets: the unbroken conservation laws $l^{\lambda_{u}}$ and the broken conservation laws $l^{\lambda_{b}}$. The unbroken conservation laws are left null eingenvectors of the submatrix of the internal species, namely $l^{\lambda_{u}}_{\alpha_{x}}\mathbb{S}^{\alpha_{x}}_{\rho}=0$ and $l^{\lambda_{u}}_{\alpha_{y}}=0$, whereas the broken conservation laws are not, namely $l^{\lambda_{b}}_{\alpha_{x}}\mathbb{S}^{\alpha_{x}}_{\rho}\neq 0$ for at least one $\rho\in\mathcal{R}$. Therefore, the quantities

[TABLE]

are conserved even if the RDS is open, $\mathrm{d}_{t}L^{\lambda_{u}}=0$, unlike $L^{\lambda_{b}}=\int_{V}\mathrm{d}\boldsymbol{r}\text{ }l^{\lambda_{b}}_{\alpha}Z^{\alpha}$. They are called unbroken conserved quantities. Notice that chemostating a species does not always break a conservation law Rao and Esposito (2016). We thus distinguish the set of controlled species $\alpha_{y_{b}}\in Y_{b}\subseteq Y$ breaking all the broken conservation laws from the others $\alpha_{y_{p}}\in Y_{p}=Y\setminus Y_{b}$. This allows us to introduce the so-called moieties

[TABLE]

where $\hat{l}^{\alpha_{y_{b}}}_{\lambda_{b}}$ denotes the elements of the inverse of the matrix whose entries are $l^{\lambda_{b}}_{\alpha_{y_{b}}}$ (see Ref. Falasco et al. (2018) for details). The moieties represent the concentration of parts of molecules which are exchanged with the environment through the chemostats since their time evolution is determined only by the external currents, $\mathrm{d}_{t}M^{\alpha_{y_{b}}}(t)=\hat{l}^{\alpha_{y_{b}}}_{\lambda_{b}}\int_{V}\mathrm{d}\boldsymbol{r}\text{ }l^{\lambda_{b}}_{\alpha_{y}}I^{\alpha_{y}}(\boldsymbol{r},t)$.

The thermodynamic equilibrium is characterized by vanishing reaction currents $j^{\rho}=0$, diffusion currents $\boldsymbol{J}^{\alpha}=0$, and external currents $I^{\alpha_{y}}=0$. The equilibrium concentrations $Z_{\text{eq}}^{\alpha}$ are, consequently, homogeneously distributed over the volume $V$. The existence of such equilibrium, together with mass-action kinetics, implies the so-called local detailed balance condition for the kinetic constants $k_{\pm\rho}$ of the chemical reactions,

[TABLE]

This constraint is assumed to be valid to ensure that closed systems relax to equilibrium and hence thermodynamic consistency Polettini and Esposito (2014). Nonequilibrium conditions can be created by chemostatting certain species. These may be homogeneously distributed or not as well as constant in time or not. Our description of the reaction-diffusion dynamics presumes that all degrees of freedom other than concentrations in space to be equilibrated, as the temperature $T$ and the pressure of the solvent. In this way, thermodynamic state functions can be specified by their equilibrium form but expressed in terms of nonequilibrium concentrations, like for the chemical potential $\mu_{\alpha}(\boldsymbol{r},t)={\mu}_{\alpha}^{\circ}+RT\ln(Z_{\alpha}(\boldsymbol{r},t))$ of each species (${\mu}_{\alpha}^{\circ}$ is the standard chemical potential). In this respect, the local detailed balance condition can be restated establishing a correspondence between the kinetic constants $k_{\pm\rho}$, namely the dynamics of the RDSs, and the standard chemical potentials $\mu_{\alpha}^{\circ}$, namely the thermodynamics:

[TABLE]

We consider now the nonequilibrium thermodynamic characterization of RDSs. First, we recapitulate the global thermodynamic theory developed in Ref. Falasco et al. (2018) and then we generalize it to a local formulation. The global second law of thermodynamics for open RDSs can be stated as the following balance equation

[TABLE]

for the so-called semigrand Gibbs free energy $\mathcal{G}(t)$. $T\dot{\Sigma}$ is the non-negative entropy production rate accounting for the energy dissipation, $\dot{W}_{\text{driv}}^{\mathcal{G}}$ is the driving work rate needed to manipulate the concentration of the chemostatted species in a time dependent way and $\dot{W}_{\text{nc}}$ is the nonconservative work rate spent to prevent the system from relaxing towards equilibrium. We examine in detail each term in Eq. (9). Whenever possible, we provide the expression of the thermodynamic quantities according to their densities which will then be used in the local formulation of thermodynamics. With a slight abuse of notation, we will use the same name for the density and the corresponding global quantities. For example, we will refer to both $T\dot{\sigma}$ and $T\dot{\Sigma}=\int_{V}\mathrm{d}\boldsymbol{r}\text{ }T\dot{\sigma}$ as entropy production rate.

In analogy to equilibrium thermodynamics when passing from canonical to grand canonical ensembles, the semigrand Gibbs free energy

[TABLE]

is obtained from the Gibbs free energy $G$ by eliminating the energetic contributions of the matter exchanged with the reservoirs. The latter amounts to the moieties $M^{\alpha_{y_{b}}}$ of Eq. (6), times the reference values of their chemical potentials $\mu^{\text{ref}}_{\alpha_{y_{b}}}$ which are the values of chemical potential fixed by the chemostats $\alpha_{y_{b}}$. If a chemostat $\alpha_{y_{b}}$ sets different values of the chemical potential $\mu_{\alpha_{y_{b}}}(\boldsymbol{r})$ for different points $\boldsymbol{r}$ in the volume $V$, then the reference chemical potential can be chosen arbitrarily among these values. This is equivalent to set the minimum value of $\mathcal{G}$ as it will become clear later. The Gibbs free energy $G(t)=\int_{V}\mathrm{d}\boldsymbol{r}\text{ }g(\boldsymbol{r},t)$ is the integral over the volume $V$ of the Gibbs free energy density $g(\boldsymbol{r},t)$ of ideal dilute solutions

[TABLE]

with $RTZ^{\mathcal{S}}=RT\sum_{\alpha}Z^{\alpha}$ accounting for the contribution of the solvent Ge and Qian (2016); Fermi (1956).

The total entropy production rate $\dot{\Sigma}$ consists of two parts $\dot{\Sigma}=\dot{\Sigma}_{\text{rct}}+\dot{\Sigma}_{\text{diff}}$, i.e., the reaction $\dot{\Sigma}_{\text{rct}}$ and the diffusion $\dot{\Sigma}_{\text{diff}}$ part:

[TABLE]

with

[TABLE]

The free energy of reaction $-(\mu_{\alpha}\mathbb{S}^{\alpha}_{\rho})$ and the variation of the chemical potential across space $-\nabla\mu_{\alpha}$ are thermodynamic forces driving the system towards equilibrium. The densities of the entropy production rates, as well as their global counterparts, are non-negative quantities and they vanish only at equilibrium. Indeed, the reaction entropy production rate can be written as $T\dot{\sigma}_{\text{rct}}=RTj^{\rho}\ln(j_{+\rho}/j_{-\rho})\geq 0$, with the aid of the local detailed balance (8) and the mass-action kinetics (2). Similarly, the diffusion entropy production rate can be written as $T\dot{\sigma}_{\text{diff}}=RTD^{\alpha}|\nabla Z_{\alpha}|^{2}/Z_{\alpha}\geq 0$. The driving work rate takes into account the time-dependent manipulation of $\mu_{\alpha_{y_{b}}}^{\text{ref}}$ performed by the chemostats $\alpha_{y_{b}}$

[TABLE]

It obviously vanishes in autonomous systems. The nonconservative work rate

[TABLE]

where the corresponding density is given by

[TABLE]

quantifies the energetic cost of sustaining fluxes of chemical species among the chemostats by means of the forces $\mathscr{F}_{\alpha_{y}}=(\mu_{\alpha_{y}}-\mu^{\text{ref}}_{\alpha_{y_{b}}}\hat{l}^{\alpha_{y_{b}}}_{\lambda_{b}}l^{\lambda_{b}}_{\alpha_{y}})$. These forces have different origin depending on whether we consider the chemostatted species $\alpha_{y_{b}}$ or $\alpha_{y_{p}}$. Indeed, $\mathscr{F}_{\alpha_{y_{b}}}(\boldsymbol{r})=\mu_{\alpha_{y_{b}}}(\boldsymbol{r})-\mu^{\text{ref}}_{\alpha_{y_{b}}}$ is the difference of chemical potential of the same species imposed in different points of the space. It vanishes in case of homogeneous chemostatting. In contrast, $\mathscr{F}_{\alpha_{y_{p}}}(\boldsymbol{r})=\mu_{\alpha_{y_{p}}}(\boldsymbol{r})-\mu^{\text{ref}}_{\alpha_{y_{b}}}\hat{l}^{\alpha_{y_{b}}}_{\lambda_{b}}l^{\lambda_{b}}_{\alpha_{y}}$ is the difference of chemical potentials imposed by different chemostats. It vanishes if all the chemostats break a conservation law, i.e., $Y_{b}=Y$.

In the absence of nonconservative forces $\mathscr{F}_{\alpha_{y}}=0$ $\forall\alpha_{y}\in Y$, no nonconservative work is performed $\dot{W}_{\mathrm{n}c}=0$ and the system is said to be detailed balanced. This occurs when all the chemostats break a conservation law $Y_{b}=Y$, and the chemostatted species are homogeneously distributed $\mu_{\alpha_{y_{b}}}(\boldsymbol{r})=\mu_{\alpha_{y_{b}}}^{\text{ref}}$ $\forall\boldsymbol{r}\in V$. If a detailed balanced system is also autonomous $\dot{W}_{\text{driv}}=0$, it relaxes towards equilibrium. This follows from the fact that the semigrand Gibbs free energy approaches its equilibrium value $\mathcal{G}_{\text{eq}}$ monotonously in time: i) the time derivative of $\mathcal{G}(t)$ is always negative,

[TABLE]

and ii) $\mathcal{G}(t)$ is lower bounded by its equilibrium value $\mathcal{G}_{\text{eq}}$. Indeed, $\mathcal{G}$ can be expressed in terms of its reference equilibrium value $\mathcal{G}_{\text{eq}}$ as

[TABLE]

by employing the relative entropy for the concentration fields (i.e., non-normalized concentration distributions)

[TABLE]

which quantifies the difference between the two concentration distributions $\{Z^{\alpha}\}$ and $\{Z^{\alpha}_{\text{eq}}\}$. Since the relative entropy is always positive, unless $Z^{\alpha}=Z^{\alpha}_{\text{eq}}$, the semigrand Gibbs free energy is greater than or equal to its equilibrium counterpart $\mathcal{G}\geq\mathcal{G}_{\text{eq}}$. In other words, $\mathcal{G}(t)$ acts as a Lyapunov function. The reference equilibrium condition

[TABLE]

is specified according to the reference chemical potentials $\mu^{\text{ref}}_{\alpha}$. We have already introduced the reference chemical potentials $\mu^{\text{ref}}_{\alpha_{y_{b}}}$ for the controlled species $\alpha_{y_{b}}$ in Eq. (10). Here, we call the chemical potentials $\mu_{\alpha_{x}}$ and $\mu_{\alpha_{y_{b}}}$ evaluated at the equilibrium condition reference chemical potentials because they depend on $\mu^{\text{ref}}_{\alpha_{y_{b}}}$. A detailed discussion in Appendix A shows that, for the $Y_{p}$ controlled species, $\mu_{\alpha_{y_{p}}}^{\text{ref}}$ is given by $\mu_{\alpha_{y_{p}}}^{\text{ref}}=\mu_{\alpha_{y_{b}}}^{\text{ref}}\hat{l}^{\alpha_{y_{b}}}_{\lambda_{b}}l^{\lambda_{b}}_{\alpha_{y}}$. This also provides an interpretation of the thermodynamic force $\mathscr{F}_{\alpha_{y_{p}}}(\boldsymbol{r})=\mu_{\alpha_{y_{p}}}(\boldsymbol{r})-\mu^{\text{ref}}_{\alpha_{y_{b}}}\hat{l}^{\alpha_{y_{b}}}_{\lambda_{b}}l^{\lambda_{b}}_{\alpha_{y}}$ as the difference between the chemical potential of the species $\alpha_{y_{b}}$ imposed by the chemostat $\alpha_{y_{b}}$ and the corresponding reference chemical potential. In contrast, $\mu_{\alpha_{x}}^{\text{ref}}$ is specified by both the reference chemical potentials $\mu_{\alpha_{y_{b}}}^{\text{ref}}$ and the unbroken conserved quantities $L^{\lambda_{u}}$ of Eq. (5) but an explicit expression cannot be obtained. If all the conservation laws are broken, there are no more conserved quantities and also $\mu_{\alpha_{x}}^{\text{ref}}$ is given by $\mu_{\alpha_{x}}^{\text{ref}}=\mu_{\alpha_{y_{b}}}^{\text{ref}}\hat{l}^{\alpha_{y_{b}}}_{\lambda_{b}}l^{\lambda_{b}}_{\alpha_{x}}$. It is import to understand that the choice of $\mu^{\text{ref}}_{\alpha_{y_{b}}}$ of Eq. (10) set the equilibrium condition $\{Z^{\text{eq}}_{\alpha}\}$ reached by the system if it was detailed balanced. Consequently, the choice of $\mu^{\text{ref}}_{\alpha_{y_{b}}}$ set $\mathcal{G}_{\text{eq}}$ that is the minimum value of $\mathcal{G}(t)$.

We now generalize this description to a local thermodynamic theory. To this aim, a fundamental step is the identification of the proper thermodynamic potential characterizing the local dynamics of RDSs as $\mathcal{G}$ does for the global dynamics. Therefore, we seek a free energy density that i) is lower bounded by its equilibrium value and ii) approaches asymptotically its equilibrium value for undriven detailed balanced systems. One might wonder whether the semigrand Gibbs free energy density

[TABLE]

such that $\mathcal{G}(t)=\int_{V}\mathrm{d}\boldsymbol{r}\text{ }\mathpzc g(\boldsymbol{r},t)$, plays this role. Its balance equation is specified as

[TABLE]

where $T\dot{\sigma}$ is the entropy production rate of Eq. (13), $\dot{w}_{\text{nc}}$ is the nonconservative work rate of Eq. (16), and the driving work rate is given by

[TABLE]

Compared to the evolution equation (9) for the semigrand Gibbs free energy, a flow term

[TABLE]

arises at the local level. This term is absent at the global level because the diffusion currents vanish at the boundaries $\partial V$ and, consequently, $\int_{V}\mathrm{d}\boldsymbol{r}\text{ }\nabla\cdot\boldsymbol{J}^{\mathpzc g}=0$ .

The semigrand Gibbs free energy density $\mathpzc g(\boldsymbol{r},t)$ is not the proper thermodynamic potential because it is not lower bounded by its equilibrium value $\mathpzc g_{\text{eq}}$: $\mathpzc g(\boldsymbol{r},t)$ can not be written as $\mathpzc g_{\text{eq}}$ plus the relative entropy of Eq. (19). Consequently, $\mathpzc g(\boldsymbol{r},t)$ is not minimized in undriven detailed balanced systems. This is a direct consequence of the definition (21) where $\mathpzc g(\boldsymbol{r},t)$ is constructed from the Gibbs free energy density $g(\boldsymbol{r},t)$ by eliminating the energetic contributions of the matter exchanged with the chemostats. Indeed, the term $\hat{l}^{\alpha_{y_{b}}}_{\lambda_{b}}l^{\lambda_{b}}_{\alpha}Z^{\alpha}(\boldsymbol{r},t)$ in Eq. (21) is the local representation of the moieties of Eq. (6). However, at the local level diffusion allows all the species to be exchanged with neighboring regions of space which play the same role as the chemostats from a local standpoint. This means that the proper thermodynamic potential must be specified as the following free energy density

[TABLE]

where the energetic contributions of all the species exchanged through diffusion are removed from the Gibbs free energy density $g(\boldsymbol{r},t)$ of Eq. (11). The local equivalent to the term $\mu^{\text{ref}}_{\alpha_{y_{b}}}M^{\alpha_{y_{b}}}$ in Eq. (10) is precisely $\mu^{\text{ref}}_{\alpha}Z^{\alpha}$. The reference chemical potentials $\mu^{\text{ref}}_{\alpha}$ have the same meaning as in Eq. (20). If all the conservation laws are broken by the chemostats, the definition of the free energy density $\mathbbm{g}(\boldsymbol{r},t)$ in Eq. (25) is equivalent to the definition of the semigrand Gibbs free energy density $\mathpzc g(\boldsymbol{r},t)$ in Eq. (21) since $\mu^{\text{ref}}_{\alpha_{x}}Z^{\alpha}=\mu^{\text{ref}}_{\alpha_{y_{b}}}\hat{l}^{\alpha_{y_{b}}}_{\lambda_{b}}l^{\lambda_{b}}_{\alpha_{x}}Z^{\alpha_{x}}$ and $\mu^{\text{ref}}_{\alpha_{y_{p}}}Z^{\alpha}=\mu^{\text{ref}}_{\alpha_{y_{b}}}\hat{l}^{\alpha_{y_{b}}}_{\lambda_{b}}l^{\lambda_{b}}_{\alpha_{y_{p}}}Z^{\alpha_{y_{p}}}$.

We can now verify that $\mathbbm{g}(\boldsymbol{r},t)$ can be expressed in terms of the reference equilibrium free energy $\mathbbm{g}_{\text{eq}}$ as

[TABLE]

by employing the relative entropy of Eq. (19). The free energy density is therefore always greater than or equal to the equilibrium counterpart $\mathbbm{g}\geq\mathbbm{g}_{\text{eq}}$. The fact that the relative entropy of information theory Cover and Thomas (2012) appears both at the global and the local level supports the idea that thermodynamics plays a fundamental role in the characterization of the information codified in complex patterns of chemical concentrations.

Furthermore, the free energy density $\mathbbm{g}(\boldsymbol{r},t)$ satisfies the balance equation

[TABLE]

which constitutes the local formulation of the second law. The non-negative entropy production $T\dot{\sigma}$ and the nonconservative work rate $\dot{w}_{\text{nc}}$ are the same as in Eq (13) and Eq. (16), respectively. However, the driving chemical work rate

[TABLE]

now takes into account the time evolution of all $\mu^{\text{ref}}_{\alpha}$ through the time-dependent manipulation of $\mu^{\text{ref}}_{\alpha_{y_{b}}}$ performed by the chemostats $\alpha_{y_{b}}$. Indeed, the reference chemical potentials $\mu^{\text{ref}}_{\alpha_{y_{p}}}$ and $\mu^{\text{ref}}_{\alpha_{x}}$ depend on $\mu^{\text{ref}}_{\alpha_{y_{b}}}$’s (see Appendix A). The flow of free energy,

[TABLE]

describes the effects on the evolution of the free energy density due to the local diffusion currents. Note that the diffusion entropy production rate $T\dot{\sigma}_{\text{diff}}$ (resp. the free energy flow $\boldsymbol{J}^{\mathbbm{g}}$) can be further split into a contribution due to the internal species $T\dot{\sigma}_{\text{diff}}^{X}=-(\nabla\mu_{\alpha_{x}})\cdot\boldsymbol{J}^{\alpha_{x}}$ (resp. $\boldsymbol{J}^{\mathbbm{g}_{X}}=(\mu_{\alpha_{x}}-\mu_{\alpha_{x}}^{\text{ref}})\boldsymbol{J}^{\alpha_{x}}$) and one due to the controlled species $T\dot{\sigma}_{\text{diff}}^{Y}=-(\nabla\mu_{\alpha_{y}})\cdot\boldsymbol{J}^{\alpha_{y}}$ (resp. $\boldsymbol{J}^{\mathbbm{g}_{Y}}=(\mu_{\alpha_{y}}-\mu_{\alpha_{y}}^{\text{ref}})\boldsymbol{J}^{\alpha_{y}}$). We will make use of this decomposition when we examine the thermodynamic description of chemical waves. Equations (25) and (27) constitute the first major result of this paper.

For undriven detailed balanced systems $\mathbbm{g}$ approaches $\mathbbm{g}_{\text{eq}}$ asymptotically. Indeed, the global free energy,

[TABLE]

is greater than its reference equilibrium value $\mathbbm{G}\geq\mathbbm{G}_{\text{eq}}$ (because $\mathbbm{g}\geq\mathbbm{g}_{\text{eq}}$) and its time derivative is always negative

[TABLE]

since $\dot{w}^{\mathbbm{g}}_{\text{driv}}=\dot{w}_{\text{nc}}=0$ and $\int_{V}\mathrm{d}\boldsymbol{r}\text{ }\nabla\cdot\boldsymbol{J}^{\mathbbm{g}}=0$. In other words, $\mathbbm{G}\to\mathbbm{G}_{\text{eq}}$ monotonously in time. While this implies that $\mathbbm{g}$ approaches $\mathbbm{g}_{\text{eq}}$ asymptotically, i.e., $\mathbbm{g}$ is minimized, it is not granted that its time derivative $\partial_{t}\mathbbm{g}=-T\dot{\sigma}-\nabla\cdot\boldsymbol{J}^{\mathbbm{g}}$ is always negative because of the free energy flow. Consider for example the case where the concentration for all the species except one is equal to its equilibrium value at one specific point $\boldsymbol{r}^{\prime}$ in space: $Z^{\alpha}(\boldsymbol{r},t)=Z^{\alpha}_{\text{eq}}$ $\forall\alpha\neq\alpha^{\prime}$ if and only if $\boldsymbol{r}=\boldsymbol{r}^{\prime}$. The concentration of the species $\alpha^{\prime}$ differs instead from its equilibrium value for an arbitrarily small positive number, $Z^{\alpha^{\prime}}(\boldsymbol{r}^{\prime},t)-Z^{\alpha^{\prime}}_{\text{eq}}=\epsilon>0$. Expanding the time derivative of the free energy $\partial_{t}\mathbbm{g}(\boldsymbol{r},t)$ for $\boldsymbol{r}=\boldsymbol{r}^{\prime}$ in powers of $\epsilon$ and truncating the expansion at the lowest order gives $\partial_{t}\mathbbm{g}(\boldsymbol{r}^{\prime},t)=\epsilon\mathrm{d}_{\epsilon}(T\dot{\sigma}_{\text{rct}})+\epsilon D^{\alpha^{\prime}}\nabla^{2}Z^{\alpha^{\prime}}/Z^{\alpha^{\prime}}_{\text{eq}}$, with the derivative of the reaction entropy production rate evaluated for $\epsilon=0$. Therefore, the free energy increases locally ($\partial_{t}\mathbbm{g}(\boldsymbol{r}^{\prime},t)>0$) if, for instance, the curvature of the concentration field is very large at $\boldsymbol{r}=\boldsymbol{r}^{\prime}$ (e.g., $\nabla^{2}Z^{\alpha^{\prime}}=Z^{\alpha^{\prime}}_{\text{eq}}/\epsilon^{2}$) since the contribution of the reaction entropy production rate becomes negligible.

III Chemical waves

We now examine the thermodynamic description of chemical waves considered here as particular pattern solutions of RDSs also called traveling waves. The concentration of a species $Z^{\alpha}(\boldsymbol{r},t)$ evolves as a traveling wave if it propagates in space with a constant velocity $\boldsymbol{c}^{\alpha}$ and without changing its shape:

[TABLE]

We label the spatial coordinate in the traveling wave $\alpha$ frame of reference as $\tilde{\boldsymbol{r}}^{\alpha}:={\boldsymbol{r}}-{\boldsymbol{c}}^{\alpha}t$.

III.1 Dynamics of Traveling Waves

We first discuss a set of conditions that the reaction-diffusion equation (1) needs to satisfy to allow for traveling wave solutions. We will then employ these conditions in Subs. III.2 to specialize the general thermodynamic description of Sec. II for traveling waves. We start by splitting the chemical species into two disjoint subsets, i.e., the traveling waves $\alpha_{x}\in X$ and the other species $\alpha_{y}\in Y$. The traveling waves are internal species by construction. We use the same notation for the nontraveling wave species as for the controlled species in Sec. II since we will show that these former must be undriven chemostatted species.

Using Eq. (32), the reaction-diffusion equation (1) for the traveling wave $\alpha_{x}$ becomes

[TABLE]

In the frame of reference of the traveling wave $\alpha_{x}$, the terms $\nabla^{2}Z^{\alpha_{x}}$ and $\boldsymbol{c}^{\alpha_{x}}\cdot\nabla Z^{\alpha_{x}}$ are time-independent and, consequently, the above equation admits a solution if the concentration variation due to the chemical reactions is time-independent as well: $\partial_{t}(\mathbb{S}^{\alpha_{x}}_{\rho}j^{\rho}(\tilde{\boldsymbol{r}}^{\alpha_{x}}+\boldsymbol{c}^{\alpha_{x}}t,t))=0$. This condition must be satisfied for every $\tilde{\boldsymbol{r}}^{\alpha_{x}}$ and every $t$ justifying the hypothesis that each reaction current $j^{\rho}$ is time-independent in the frame of reference of the traveling wave $\alpha_{x}$, i.e., $\partial_{t}(j^{\rho}(\tilde{\boldsymbol{r}}^{\alpha_{x}}+\boldsymbol{c}^{\alpha_{x}}t,t))=0$ $\forall\rho\in\mathcal{R}$. In other words, each reaction current evolves like the traveling wave $\alpha_{x}$ with the propagation velocity $\boldsymbol{c}^{\alpha_{x}}$. Since this must be true for every traveling species, the reaction currents $j^{\rho}$ have to be time-independent in every frame of reference $\tilde{\boldsymbol{r}}^{\alpha_{x}}$. Therefore, we conclude that all the traveling waves evolve with the same velocity $\boldsymbol{c}$ and the corresponding comoving coordinate is $\tilde{\boldsymbol{r}}=\boldsymbol{r}-\boldsymbol{c}t$.

This implies that the condition $\partial_{t}(j^{\rho}(\tilde{\boldsymbol{r}}+\boldsymbol{c}t,t))=0$, using the mass-action kinetics Eq. (2) and $\partial_{t}(Z^{\alpha_{x}}(\tilde{\boldsymbol{r}}))=0$, becomes

[TABLE]

where $\nabla Z^{\alpha_{y}}$ and $\partial_{t}Z^{\alpha_{y}}$ are evaluated in frame of reference $\tilde{\boldsymbol{r}}$. The above equation must hold for every reaction current $j^{\rho}$ of every point $\tilde{\boldsymbol{r}}$ and of every time $t$. We thus assume that each term between square bracket vanish independently, $\boldsymbol{c}\cdot\nabla Z^{\alpha_{y}}+\partial_{t}Z^{\alpha_{y}}=0$. This means that the species $Y$ are either traveling waves $\partial_{t}Z^{\alpha_{y}}=-\boldsymbol{c}\cdot\nabla Z^{\alpha_{y}}$ or undriven controlled species $\partial_{t}Z^{\alpha_{y}}=0$ with a concentration field that is invariant under space translation along the propagation direction $\boldsymbol{c}\cdot\nabla Z^{\alpha_{y}}=0$. Since the species $Y$ are not traveling waves by hypothesis, they are undriven chemostatted species. In this respect, their reaction-diffusion equation simplifies to

[TABLE]

with vanishing diffusion currents in the propagation directions $\boldsymbol{c}\cdot\boldsymbol{J}^{\alpha_{y}}=-D^{\alpha_{y}}\boldsymbol{c}\cdot\nabla Z^{\alpha_{y}}=0$.

In conclusion, to allow for traveling waves to possibly arise, we consider only RDSs that satisfy the following requirements:

• (i)

the chemical species are split into two disjoint subsets, i.e., the traveling waves $X$ (internal species) and the controlled species $Y$;

• (ii)

the traveling waves are characterized by the same propagation velocity $\boldsymbol{c}$ and they have the same comoving coordinate $\tilde{\boldsymbol{r}}=\boldsymbol{r}-\boldsymbol{c}t$;

• (iii)

the concentration field of the controlled species $Y$ is invariant under space translation along the propagation direction $\boldsymbol{c}\cdot\nabla Z^{\alpha_{y}}=0$ and it is kept constant by the chemostats $\partial_{t}Z^{\alpha_{y}}=0$ (no driving work is performed, $\dot{w}_{\text{driv}}=0$).

Although we do not expect it, we note that rigorously speaking weaker conditions for traveling wave solutions may exist as our argumentation is not a mathematical proof.

III.2 Thermodynamics of Traveling Waves

We now examine how the thermodynamic quantities such as the local free energy, the entropy production rate, and the nonconservative work rate evolve when traveling waves propagate in RDSs.

Using Eq. (32), the time evolution of the free energy is given by

[TABLE]

where one can recognize a drift contribution for every traveling wave $\boldsymbol{c}\cdot\nabla Z^{\alpha_{x}}$. In other words, the free energy is “dragged” by the propagating traveling waves in the $\boldsymbol{c}$ direction.

The nonconservative work rate and the entropy production rate account for the energetic cost of supporting the wave propagation. Using Eq. (16) and Eq. (35), $\dot{w}_{\text{nc}}$ becomes

[TABLE]

The first term on the r.h.s., $-(\mu_{\alpha_{y}}-\mu^{\text{ref}}_{\alpha_{y}})\mathbb{S}^{\alpha_{y}}_{\rho}j^{\rho}$, is the energetic cost for balancing with the chemostats the concentration variations of the $Y$ species due to the chemical reactions. On the other hand, the contribution to the free energy flow $\boldsymbol{J}^{\mathbbm{g}_{Y}}$ and to the diffusion entropy production rate $T\dot{\sigma}_{\text{diff}}^{Y}$ resulting from the controlled species $\alpha_{y}$ take into account the amount of work needed to prevent modifications to the pattern profiles $Z^{\alpha_{y}}(\boldsymbol{r})$ because of diffusion.

Finally, with the aid of Eqs. (27), (36) and (37), we find that the entropy production rate reads

[TABLE]

This general result highlights explicitly which phenomena are responsible for the dissipation: i) a drift contribution due to the propagation of the waves $(\mu_{\alpha_{x}}-\mu^{\text{ref}}_{\alpha_{x}})\boldsymbol{c}\cdot\nabla Z^{\alpha_{x}}=-\partial_{t}\mathbbm{g}$, ii) the free energy flow due to the wave species $\boldsymbol{J}^{\mathbbm{g}_{X}}=(\mu_{\alpha_{x}}-\mu_{\alpha_{x}}^{\text{ref}})\boldsymbol{J}^{\alpha_{x}}$, iii) the reaction consumption/production of the chemostatted species sustaining the wave dynamics $(\mu_{\alpha_{y}}-\mu^{\text{ref}}_{\alpha_{y}})\mathbb{S}^{\alpha_{y}}_{\rho}j^{\rho}$ and iv) the diffusion of the controlled species $T\dot{\sigma}_{\text{diff}}^{Y}$. The first two contributions take into account the dynamics of the waves, while the others represent the dissipation due to the chemostatting of the controlled species. Equation (38) constitutes the second major result of this paper.

The above expressions of the thermodynamic quantities explicitly account for every process occurring while waves propagate. Furthermore, depending on the physical properties of the system of interest, they can be effectively simplified. For example, if the controlled species were homogeneously distributed $Z^{\sigma_{y}}(\boldsymbol{r})=\overline{Z}^{\sigma_{y}}$ $\forall\boldsymbol{r}\in V$, no work is required to prevent changes in the controlled species concentrations because of the diffusion. In this case, the nonconservative work rate simplifies to $\dot{w}_{\text{nc}}=-(\mu_{\alpha_{y}}-\mu^{\text{ref}}_{\alpha_{y}})\mathbb{S}^{\alpha_{y}}_{\rho}j^{\rho}$, while the diffusion entropy production rate $T\dot{\sigma}^{Y}_{\text{diff}}$ vanishes. In case the system is detailed balanced, no nonconservative forces are applied, namely $\mu_{\alpha_{y}}=\mu^{\text{ref}}_{\alpha_{y}}$, and the dissipation is only due to the propagation of the waves: $T\dot{\sigma}$ depends on the drift and the flow contribution, $T\dot{\sigma}=(\mu_{\alpha_{x}}-\mu^{\text{ref}}_{\alpha_{x}})\boldsymbol{c}\cdot\nabla Z^{\alpha_{x}}-\nabla\cdot J^{\mathbbm{g}_{X}}$. In this case, the global dissipation is solely caused by the transport mechanism of the wave dynamics $T\dot{\Sigma}=\int_{V}\mathrm{d}\boldsymbol{r}\text{ }(\mu_{\alpha_{x}}-\mu^{\text{ref}}_{\alpha_{x}})\boldsymbol{c}\cdot\nabla Z^{\alpha_{x}}$.

IV Model Systems

We apply here the thermodynamic description developed in the previous section to study the properties of the traveling waves emerging in two model systems. The first model, the Fisher-Kolmogorov equation, is detailed balanced, and the relaxation towards equilibrium occurs via a step-like wave. On the other hand, for the Brusselator model two oscillating traveling waves are maintained by nonconservative forces.

IV.1 Fisher-Kolmogorov equation

A minimal RDS displaying traveling wave solutions is the so-called Fisher-Kolmogorov equation in one infinite spatial dimension Murray (2002). It can be used to describe an autocatalytic reaction between an internal species \chX and a controlled species \chY according to the chemical equation \chX + Y ¡=¿[ $k_{+1}$ ][ $k_{-1}$ ] 2 X. The system is detailed balanced (no nonconservative forces) since one single species is chemostatted. The reaction-diffusion equation

[TABLE]

with $y$ the homogeneous concentration of the controlled species and $D$ the diffusion coefficient of the internal species, specifies the dynamics of the concentration $Z^{\alpha_{x}}(r,t)=x(r,t)$ of the internal species. Waves solutions $x(r,t)=x(r-ct)$ emerge for propagation velocities $c$ greater than or equal to a critical value $c\geq\overline{c}=2\sqrt{Dk_{+1}y}$, and they are all characterized by a step-like profile which formally means i) $\lim_{\tilde{r}\to+\infty}x(\tilde{r})=0$, ii) $\lim_{\tilde{r}\to-\infty}x(\tilde{r})=x_{\text{eq}}$ and iii) $\partial_{\tilde{r}}x<0$ $\forall\tilde{r}$. Note that $\tilde{r}$ denotes the single spatial coordinate in the frame of reference of the traveling wave, $\tilde{r}=r-ct$, while $x_{\text{eq}}$ labels the equilibrium concentration of the internal species, $x_{\text{eq}}=k_{+1}y/k_{-1}$, corresponding to the amplitude of the wave. According to these general properties, the thermodynamic quantities of the global system, namely the entropy production rate $T\dot{\Sigma}$ and the time derivative of the free energy $\mathrm{d}_{t}\mathbbm{G}$ (obtained by integrating the expressions in Eq. (36) and (38) over the volume $V$), can be directly related to the features of the traveling wave, namely the amplitude $x_{\text{eq}}$ and the propagation velocity $c$. Indeed, by considering the infinite volume limit, one obtains

[TABLE]

establishing a linear dependence between the thermodynamic quantities and the dynamical ones. The above equation, valid for every traveling wave solution, proves that the higher the amplitude and/or the velocity, the greater the dissipation of the global system during the propagation.

Particular traveling wave solutions have to be taken into account if one wants to investigate the local profile of the thermodynamic quantities. For the Fisher-Kolmogoroff equation, both asymptotic solutions and one particular analytical solution are available Murray (2002). Since the traveling waves share the same general behavior, we consider the analytical solution

[TABLE]

with propagation velocity $c=(5/\sqrt{6})\sqrt{Dk_{+1}y}$, to calculate the entropy production rate $T\dot{\sigma}$, the free energy flow $\partial_{r}J^{\mathbbm{g}}$ and the time derivative of the free energy $\partial_{t}\mathbbm{g}$ according to their definitions given in the Eqs. (38), (29) and (36).

By inspecting Fig. 1 where these thermodynamic quantities are displayed, two main features can be noticed. First of all, the variation of the thermodynamic quantities occurs around the wavefront of the concentration field as one would expect. Secondly, the entropy production rate is larger than the absolute value of the free energy flow, $T\dot{\sigma}>|\partial_{r}J^{\mathbbm{g}}|$, granting that the free energy continuously approaches its equilibrium value even locally, i.e., $\partial_{t}\mathbbm{g}<0$ $\forall\tilde{r}$. In other words, the system is relaxing continuously towards equilibrium for every point of the traveling wave profile. In this model system, the propagation dynamics of the traveling wave does not require to be sustained by nonconservative forces and, consequently, it results to be particularly effective from this specific thermodynamic point of view.

IV.2 Brusselator model

A simple nondetailed balanced system where oscillating traveling waves emerge is the Brusselator model in one spatial dimension Prigogine and Lefever (1968). It describes the transformation of two internal species, an activator \chX_1 and an inhibitor \chX_2, according to the chemical equations

[TABLE]

which involve four chemostatted species \chY_1, \chY_2, \chY_3, \chY_4. We choose \chY_1 and \chY_2 as the reference chemostatted species breaking the two conservation laws and, consequently, determining the equilibrium condition Falasco et al. (2018). The reaction-diffusion equation for the evolution in time and space $r\in V=[0,L]$ of the concentrations $Z^{\alpha_{x_{1}}}(r,t)=x_{1}(r,t)$ and $Z^{\alpha_{x_{2}}}(r,t)=x_{2}(r,t)$ of the two internal species is specified as

[TABLE]

with $y_{1},\dots,y_{4}$ the homogeneous concentrations of the controlled species, and $D_{1}$ and $D_{2}$ the diffusion coefficients of the internal species. Equation (42) admits a uniform steady-state solution $(x_{1}^{\text{ss}},x_{2}^{\text{ss}})$ that becomes unstable if the concentration $y_{2}$ of the controlled species \chY_2, used here as a bifurcation parameter, exceeds a critical value, $y_{2}>\overline{y}_{2}$. Depending on the critical point, different phenomena emerge including chemical oscillations, Turing patterns, and oscillating traveling waves Auchmuty and Nicolis (1975); Herschkowitz-Kaufman (1975); Auchmuty and Nicolis (1976). The bifurcation point leading to oscillating traveling waves with velocity $c$ and wavenumber $\nu=2\pi n/L$ (with $n\in\mathbb{N}$ the number of oscillations in the volume $V$ with periodic boundary conditions) is determined by the concentration $\overline{y}_{2}$ such that the matrix representing the linearized dynamical system (42) around the steady-state acquires the pure imaginary eigenvalues $\pm\imath\nu\overline{c}$ (with $\overline{c}$ the critical velocity).

In Ref. Auchmuty and Nicolis (1976), J. F. G. Auchmuty and G. Nicolis employed a perturbation expansion near the onset of instability, $y_{2}\simeq\overline{y}_{2}$, for identifying the oscillating wave solutions for the concentration $x_{1}(r,t)$ and $x_{2}(r,t)$ around their steady-state values in the case of irreversible chemical reactions. We use here the same analysis (summarized in Appendix B) for the reversible model by assuming that the backward reaction currents are almost negligible $j^{-\rho}\simeq 0$ $\forall\rho\in\mathcal{R}$. This means that the kinetic constants of the forward reactions are greater than the ones of the backward reactions $k_{+\rho}\gg k_{-\rho}$, and the concentrations $y_{3}$ and $y_{4}$ are significantly smaller than their equilibrium values, $y_{3,\text{eq}}=k_{+2}k_{+3}y_{2}/k_{-2}k_{-3}$ and $y_{4,\text{eq}}=k_{+1}k_{+4}y_{1}/k_{-1}k_{-4}$. In this case, the existence of traveling wave solutions of Eq. (42) with wavenumber $\nu=2\pi n/L$ is constrained by the following condition

[TABLE]

granting the existence of the propagation velocity $c$. The corresponding critical concentration is specified as

[TABLE]

Equation (43) and (44) mean that there is only a finite set of wavenumbers $\nu$ and corresponding critical concentrations $\overline{y}_{2}(\nu)$ such that the steady state becomes unstable with respect to traveling wave solutions.

Similarly to what has been done for the Fisher-Kolmogorov equation discussed in Subs. IV.1, we establish a correspondence between a global thermodynamic quantity of the Brusselator model and the features of the traveling waves without considering specific solutions. In particular, the global nonconservative work rate $\dot{W}_{\text{nc}}$, i.e., the integral over the volume $V$ of the work rate of Eq. (37), $\dot{W}_{\text{nc}}=-\int_{0}^{L}\mathrm{d}r\text{ }(\mu_{\alpha_{y}}-\mu^{\text{ref}}_{\alpha_{y}})\mathbb{S}^{\alpha_{y}}_{\rho}j^{\rho}$, can be expressed as a function of the wave wavenumber $\nu$. Indeed, by taking into account that i) the chemical potentials of the species \chY_1 and \chY_2 correspond to their reference potentials, $\mu_{\alpha_{y_{1}}}=\mu^{\text{ref}}_{\alpha_{y_{1}}}$ and $\mu_{\alpha_{y_{2}}}=\mu^{\text{ref}}_{\alpha_{y_{2}}}$, ii) the reaction currents for $\rho=3,\text{ }4$ can be specified as $j^{3}\simeq k_{+2}x_{1}y_{2}$ and $j^{4}\simeq k_{+4}x_{1}$, iii) the oscillation of each traveling wave is centered around the corresponding steady-state concentration $(x_{1}^{\text{ss}},x_{2}^{\text{ss}})\simeq((k_{+1}y_{1}/k_{+4}),(k_{+2}y_{2}/k_{+3}x_{1}^{\text{ss}}))$, the nonconservative work rate reads

[TABLE]

The concentration $y_{2}$ can be approximated with its critical value $\overline{y}_{2}(\nu)$ defined in Eq. (44) if wave solutions near the onset of the instability are considered. Since $y_{3,\text{eq}}\gg y_{3}$, the global nonconservative work rate is a monotonically increasing function of the wavenumber:

[TABLE]

This means that the energetic cost of sustaining oscillating traveling waves for the Brusselator model increases with the wavenumber $\nu$.

We consider now the specific traveling wave solutions,

[TABLE]

whose derivation is discussed in Appendix B, to study the local profile of the thermodynamic quantities defined in Eqs. (36), (37), (38) and (29), and displayed in Fig. 2 for the case with $n=1$. The particular values of the physical quantities employed for the plot, like the concentration of the chemostatted species, are reported in Appendix B.

First of all, one can verify that all the thermodynamic quantities share the same period as the traveling waves. Secondly, there is a significant difference between the nonconservative work rate and the entropy production rate, on the one hand, and the time derivative of the free energy and the divergence of its flow, on the other hand. The nonconservative work spent to sustain the propagation of the waves is mainly dissipated, $\dot{w}_{\text{nc}}\sim T\dot{\sigma}$, resulting in a limited variation of the free energy in time $\dot{w}_{\text{nc}}\sim T\dot{\sigma}\gg|\partial_{t}\mathbbm{g}|$. This behavior is caused by the need of keeping the system far from equilibrium to allow the emergence of oscillating traveling waves. From this thermodynamic point of view, the propagation of oscillating traveling waves in the Brusselator model is energetically very expensive: it requires that a lot of work is performed.

A similar general behavior is observed if one considers traveling waves emerging under different conditions, such as different concentrations of the controlled species and different wavenumber. In this regard, an interesting difference is the increase of the global nonconservative work rate with the wavenumbers according to Eq. (45). A similar trend can also be observed for the amplitude of the oscillating nonconservative work rate, as well as the amplitude of the free energy.

V Conclusions

In this work, we established a local thermodynamic theory of reaction-diffusion systems valid arbitrarily far from equilibrium. We then used it to investigate the propagation of chemical waves and its energetic cost. To this aim, the identification of the proper thermodynamic potential as well as establishing its balance equations have been the fundamental steps. We showed a clear correspondence between the thermodynamic quantities and the dynamical processes occurring while chemical waves propagate, both in general and for two model systems. The global entropy production rate of the Fisher- Kolmogorov waves is linearly dependent on the velocity and the wave amplitude. The global nonconservative work rate is a monotonically increasing function of the wavenumber for the periodic waves of the Brusselator model. The relationship between the features of the concentration profiles of the chemical waves and the thermodynamic quantities has also been analyzed in these models at the local level.

We focused here on the propagation dynamics of chemical waves, but our thermodynamic description could be used to study any other reaction-diffusion processes. A direct application might be the characterization of the energetic cost of creating wave patterns from generic initial conditions and of the response to perturbations in model systems like those studied in this work. Furthermore, we think that our approach represents the proper framework to study the correspondence between the semantic information codified in complex chemical patterns like chemical waves and the syntactic information. The semantic information represents the information content relevant for a specific system which depends on how the patterns are interpreted by for instance a receptor while the syntactic information represents the whole physical information content Kolchinsky and Wolpert (2018). These topics are becoming more important nowadays since they play a crucial role in the development of new ways of transferring and computing information with chemical reactions Rossi et al. (2008); Padirac et al. (2013); Tayar et al. (2015); Adamatzky et al. (2005). They should also shed light on how biosystems perform these operations with high efficiency Sartori et al. (2014); Barato et al. (2014); Ito and Sagawa (2015).

Finally, our work is based on a deterministic description of the reaction-diffusion dynamics, which emerges in the large size limit of a microscopic description in terms of stochastically reacting and diffusing chemicals. A thermodynamic at that level is necessary to describe intracellular signaling because of the limited number of molecules. For instance, the importance of stochasticity in intracellular calcium signaling is still debated Falcke (2003); Skupin et al. (2010).

VI Acknowledgments

This research was funded by the European Research Council project NanoThermo (ERC-2015-CoG Agreement No. 681456).

Appendix A Reference Chemical Potentials

We examine here the relations between the reference chemical potentials of different species $\mu^{\text{ref}}_{\alpha}$ for open RDSs. Note that the classification of the species $\mathcal{S}$ as internal species $X$ or controlled species breaking the broken conservation laws $Y_{b}$ or other controlled species $Y_{p}$ is particularly important in this framework. As already mentioned in the main text, the reference chemical potentials are the values of the chemical potentials at the equilibrium that would be reached if the system were detailed balanced. For this reason and the local detailed balanced condition $k_{+\rho}/k_{-\rho}=\prod_{\alpha}(Z^{\alpha}_{\text{eq}})^{\mathbb{S}^{\alpha}_{\rho}}=\exp(-\mathbb{S}^{\alpha}_{\rho}\mu_{\alpha}^{\circ}/RT)$, the set of reference chemical potentials is a left null vector of the stoichiometric matrix,

[TABLE]

This means that the reference chemical potentials $\mu^{\text{ref}}_{\alpha}$ can be written as a linear combination of the closed system conservation laws $l^{\lambda}_{\alpha}$,

[TABLE]

where $\psi_{\lambda}$ are real coefficients. Because of Eq. (49), the number of independent reference chemical potentials is equal to the number of conservation laws $|\lambda|$ of the closed RDS. The other chemical potentials are determined according to the constraints of Eq. (48).

First of all, we show that the reference chemical potentials of the chemostatted species which break all the broken conservation laws $\mu^{\text{ref}}_{\alpha_{y_{b}}}$ set the value of the reference chemical potentials of the other controlled species $\mu^{\text{ref}}_{\alpha_{y_{p}}}$. By definition of unbroken conservation laws, i.e., $l^{\lambda_{u}}_{\alpha_{x}}\mathbb{S}^{\alpha_{x}}_{\rho}=0$ $\forall\rho\in\mathcal{R}$ and $l^{\lambda_{u}}_{\alpha_{y}}=0$, the reference chemical potentials $\mu^{\text{ref}}_{\alpha_{y_{b}}}$ (resp. $\mu^{\text{ref}}_{\alpha_{y_{p}}}$) are defined as a linear combination of only the broken conservation laws, $\mu^{\text{ref}}_{\alpha_{y_{b}}}=\psi_{\lambda_{b}}l^{\lambda_{b}}_{\alpha_{y_{b}}}$ (resp. $\mu^{\text{ref}}_{\alpha_{y_{p}}}=\psi_{\lambda_{b}}l^{\lambda_{b}}_{\alpha_{y_{p}}}$). The matrix whose entries are $l^{\lambda_{b}}_{\alpha_{y_{b}}}$ is square and nonsingular so that it can be inverted. The elements of the inverse matrix are denoted by $\hat{l}_{\lambda_{b}}^{\alpha_{y_{b}}}$, and the reference chemical potentials $\mu^{\text{ref}}_{\alpha_{y_{p}}}$ can now be specified as

[TABLE]

since $\psi_{\lambda_{b}}=\mu^{\text{ref}}_{\alpha_{y_{b}}}\hat{l}_{\lambda_{b}}^{\alpha_{y_{b}}}$.

Secondly, we show that the value of reference chemical potentials for the internal species $\mu^{\text{ref}}_{\alpha_{x}}$ depends on the reference chemical potentials $\mu^{\text{ref}}_{\alpha_{y_{b}}}$ and the unbroken conserved quantities $L^{\lambda_{u}}=\int_{V}\mathrm{d}\boldsymbol{r}\text{ }l^{\lambda_{u}}_{\alpha_{x}}Z^{\alpha_{x}}$. Because of $\psi_{\lambda_{b}}=\mu^{\text{ref}}_{\alpha_{y_{b}}}\hat{l}_{\lambda_{b}}^{\alpha_{y_{b}}}$,

[TABLE]

with the coefficients $\psi_{\lambda_{u}}$ to be determined by employing the unbroken conserved quantities evaluated for the equilibrium concentrations $L^{\lambda_{u}}=Vl^{\lambda_{u}}_{\alpha_{x}}Z^{\alpha_{x}}_{\text{eq}}$. Indeed, the conserved quantities can be written as

[TABLE]

establishing $|\lambda_{u}|$ nonlinear constrains for the $|\lambda_{u}|$ coefficients $\psi_{\lambda_{u}}$. By solving the above system of equations for the coefficients $\psi_{\lambda_{u}}$, one obtains $\psi_{\lambda_{u}}$ as a function of $\{L^{\lambda_{u}}\}$ and $\{\mu^{\text{ref}}_{\alpha_{y_{b}}}\}$. Unfortunately, explicit expressions cannot be obtained from Eq. (52). If there are no conserved quantities, $\mu_{\alpha_{x}}^{\text{ref}}=\mu_{\alpha_{y_{b}}}^{\text{ref}}\hat{l}^{\alpha_{y_{b}}}_{\lambda_{b}}l^{\lambda_{b}}_{\alpha_{x}}$. In conclusion, the reference chemical potentials $\mu^{\text{ref}}_{\alpha_{y_{b}}}$ and the conserved quantities $L^{\lambda_{u}}$ set the values of all the other reference chemical potentials $\mu^{\text{ref}}_{\alpha_{x}}$, $\mu^{\text{ref}}_{\alpha_{y_{p}}}$.

Appendix B Traveling Waves in the Brusselator Model

We summarize here the derivation of oscillating traveling wave solutions of the Brusselator model in one spatial dimension $r\in[0,L]$ with irreversible chemical reactions proposed by J. F. G. Auchmuty and G. Nicolis Auchmuty and Nicolis (1976). Since the backward reaction currents are neglected, the reaction-diffusion equation (42) reads, in vector form,

[TABLE]

Oscillating traveling wave solutions of Eq. (53) are periodic functions of the only comoving coordinate: $x_{1}(r,t)=x_{1}(\tilde{r})=x_{1}(\tilde{r}\pm L)$ and $x_{2}(r,t)=x_{2}(\tilde{r})=x_{2}(\tilde{r}\pm L)$, with $\tilde{r}=r-ct$. In this regard, the steady-state $(x_{1}^{\text{ss}},x_{2}^{\text{ss}})=((k_{+1}y_{1}/k_{+4}),(k_{+2}y_{2}/k_{+3}x_{1}^{\text{ss}}))$ can be interpreted as a trivial traveling wave with null velocity. By considering the linearization of the dynamical system (53) around the steady-state according to the operator

[TABLE]

one can identify the conditions where new wave solutions may bifurcate from the steady-state. Any such a solution, written as a perturbation around the steady-state $(x_{1}^{\text{p}},x_{2}^{\text{p}})=(x_{1}-x_{1}^{\text{ss}},x_{2}-x_{2}^{\text{ss}})$, must have the following form in the linear regime

[TABLE]

with $(a_{1},a_{2})$ still unknown coefficients and $\nu=2\pi n/L$ ($n\in\mathbb{N}$), because of the periodic boundary conditions. Consequently, $(x_{1}^{\text{p}},x_{2}^{\text{p}})$ of Eq. (55) must be an eigenvector of the linear operator $\hat{\mathcal{J}}$ corresponding to the complex eigenvalue $-\imath\nu c$. To guarantee the existence of such a solution, the matrix representation of the operator $\hat{\mathcal{J}}$ on the basis element $\exp(\imath\nu\tilde{r})$

[TABLE]

must have purely complex eigenvalues meaning that its trace vanishes $\mathcal{T}=0$, while its determinant is positive $\mathcal{D}>0$. The condition $\mathcal{T}=0$ determines the critical value $\overline{y}_{2}$ of the concentration of the chemostatted species \chY_2, considered here as the bifurcation parameter, which is specified in Eq. (44). On the other hand, the constraint $\mathcal{D}>0$ is equivalent to Eq. (43) when the critical concentration $\overline{y}_{2}$ is employed, and it ensures that the value of the critical velocity $\overline{c}(\nu)$, namely the wave velocity evaluated for $y_{2}=\overline{y}_{2}$, is a real number:

[TABLE]

When the concentration $y_{2}$ exceeds its critical value $\overline{y}_{2}(\nu)$, the steady-state becomes unstable compared to the traveling waves with wavenumber $\nu$. Waves solutions of the reaction-diffusion equation close to the bifurcation point can be identified by expanding the dynamical system (53) in powers of a small parameter $\epsilon$ representing the distance from the threshold. To this aim, Eq. (53) can be written in such a way that i) the difference between the actual concentration $y_{2}$ (resp. velocity $c$) and the critical one $\overline{y}_{2}(\nu)$ (resp. $\overline{c}(\nu)$) is made explicit, and ii) the linear contributions to the dynamics are gathered together and separated from the nonlinear one. In this way, Eq. (53) reads

[TABLE]

where $\hat{\overline{\mathcal{J}}}$ is the operator $\hat{{\mathcal{J}}}$ evaluated at the critical concentration $\overline{y}_{2}(\nu)$, while $h(x_{1},x_{2})$ takes into account the nonlinear dynamics. We shall now expand every term in Eq. (58) in powers of $\epsilon$

[TABLE]

and then iteratively solve each equation obtained by gathering equal powers of $\epsilon$. Following this perturbation procedure, one derives the traveling wave solutions with wavenumber $\nu$ of the Brusselator model. The details of the perturbation expansion in power of $\epsilon$ are discussed in Ref. Auchmuty and Nicolis (1976), and they are not reported here.

The particular traveling wave solution in Eq. (47), as well as any other quantities mentioned in the following, is specified according to the set of unit measures reported in the Table 1.

It has been obtained under the hypothesis of kinetic constants $k_{+\rho}=1$ and stopping the perturbation expansion at the second order. The explicit expressions of the coefficients $\epsilon$, $\phi$, $\phi_{1}$, $\phi_{2}$, $\theta$, $\theta_{1}$, $\theta_{2}$ can be found in Ref. Auchmuty and Nicolis (1976). For the plot in Fig. 2, we adopt the following numerical values of the physical quantities: $y_{1}=0.5$, $y_{2}=1.01\overline{y}_{2}$, $D_{2}=0.05$, $y_{3}=10$, $y_{4}=1$ and $k_{-\rho}=10^{-4}$. Similar results to those discussed in the Subs. IV.2 are obtained with different numerical values as long as the concentrations of the chemostatted species \chY_3 and \chY_4 are much smaller than their equilibrium ones.

Finally, it has to be mentioned that our thermodynamic theory is based on the implicit assumption that a reference equilibrium condition is well defined. This means that it can be applied only to the Brusselator model with reversible chemical reactions. On the other hand, we employ the analytical traveling wave solutions of the irreversible Brusselator as an approximate expression of the traveling wave solutions for the case of negligible, but not null, backward reaction currents. This choice appears to be reasonable since the relations between the thermodynamic quantities, like $\partial_{t}\mathbbm{g}=-T\dot{\sigma}+\dot{w}_{\text{nc}}-\partial_{r}J^{\mathbbm{g}}$, are still verified.