A Unification between Dynamical System Theory and Thermodynamics Involving an Energy, Mass, and Entropy State Space Formalism

Abstract

In this paper, we combine the two universalisms of thermodynamics and dynamical systems theory to develop a dynamical system formalism for classical thermodynamics. Specifically, using a compartmental dynamical system energy flow model involving heat flow, work energy, and chemical reactions, we develop a state-space dynamical system model that captures the key aspects of thermodynamics, including its fundamental laws. In addition, we show that our thermodynamically consistent dynamical system model is globally semistable with system states converging to a state of temperature equipartition. Furthermore, in the presence of chemical reactions, we use the law of mass-action and the notion of chemical potential to show that the dynamic system states converge to a state of temperature equipartition and zero affinity corresponding to a state of chemical equilibrium.

1. Introduction

Thermodynamics is a physical branch of science that governs the thermal behavior of dynamical systems from those as simple as refrigerators to those as complex as our expanding universe. The laws of thermodynamics involving conservation of energy and nonconservation of entropy are, without a doubt, two of the most useful and general laws in all sciences. The first law of thermodynamics, according to which energy cannot be created or destroyed but is merely transformed from one form to another, and the second law of thermodynamics, according to which the usable energy in an adiabatically isolated dynamical system is always diminishing in spite of the fact that energy is conserved, have had an impact far beyond science and engineering. The second law of thermodynamics is intimately connected to the irreversibility of dynamical processes. In particular, the second law asserts that a dynamical system undergoing a transformation from one state to another cannot be restored to its original state and at the same time restore its environment to its original condition. That is, the status quo cannot be restored everywhere. This gives rise to a monotonically increasing quantity known as entropy. Entropy permeates the whole of nature, and unlike energy, which describes the state of a dynamical system, entropy is a measure of change in the status quo of a dynamical system.

There is no doubt that thermodynamics is a theory of universal proportions whose laws reign supreme among the laws of nature and are capable of addressing some of science’s most intriguing questions about the origins and fabric of our universe. The laws of thermodynamics are among the most firmly established laws of nature and play a critical role in the understanding of our expanding universe. In addition, thermodynamics forms the underpinning of several fundamental life science and engineering disciplines, including biological systems, physiological systems, chemical reaction systems, ecological systems, information systems, and network systems, to cite but a few examples. While from its inception its speculations about the universe have been grandiose, its mathematical foundation has been amazingly obscure and imprecise [1,2,3,4]. This is largely due to the fact that classical thermodynamics is a physical theory concerned mainly with equilibrium states and does not possess equations of motion. The absence of a state space formalism in classical thermodynamics, and physics in general, is quite disturbing and in our view largely responsible for the monomeric state of classical thermodynamics.

In recent research [4,5,6], we combined the two universalisms of thermodynamics and dynamical systems theory under a single umbrella to develop a dynamical system formalism for classical thermodynamics so as to harmonize it with classical mechanics. While it seems impossible to reduce thermodynamics to a mechanistic world picture due to microscopic reversibility and Poincaré recurrence, the system thermodynamic formulation of [4] provides a harmonization of classical thermodynamics with classical mechanics. In particular, our dynamical system formalism captures all of the key aspects of thermodynamics, including its fundamental laws, while providing a mathematically rigorous formulation for thermodynamical systems out of equilibrium by unifying the theory of heat transfer with that of classical thermodynamics. In addition, the concept of entropy for a nonequilibrium state of a dynamical process is defined, and its global existence and uniqueness is established. This state space formalism of thermodynamics shows that the behavior of heat, as described by the conservation equations of thermal transport and as described by classical thermodynamics, can be derived from the same basic principles and is part of the same scientific discipline.

Connections between irreversibility, the second law of thermodynamics, and the entropic arrow of time are also established in [4,6]. Specifically, we show a state irrecoverability and, hence, a state irreversibility nature of thermodynamics. State irreversibility reflects time-reversal non-invariance, wherein time-reversal is not meant literally; that is, we consider dynamical systems whose trajectory reversal is or is not allowed and not a reversal of time itself. In addition, we show that for every nonequilibrium system state and corresponding system trajectory of our thermodynamically consistent dynamical system, there does not exist a state such that the corresponding system trajectory completely recovers the initial system state of the dynamical system and at the same time restores the energy supplied by the environment back to its original condition. This, along with the existence of a global strictly increasing entropy function on every nontrivial system trajectory, establishes the existence of a completely ordered time set having a topological structure involving a closed set homeomorphic to the real line, thus giving a clear time-reversal asymmetry characterization of thermodynamics and establishing an emergence of the direction of time flow.

In this paper, we reformulate and extend some of the results of [4]. In particular, unlike the framework in [4] wherein we establish the existence and uniqueness of a global entropy function of a specific form for our thermodynamically consistent system model, in this paper we assume the existence of a continuously differentiable, strictly concave function that leads to an entropy inequality that can be identified with the second law of thermodynamics as a statement about entropy increase. We then turn our attention to stability and convergence. Specifically, using Lyapunov stability theory and the Krasovskii–LaSalle invariance principle [7], we show that for an adiabatically isolated system, the proposed interconnected dynamical system model is Lyapunov stable with convergent trajectories to equilibrium states where the temperatures of all subsystems are equal. Finally, we present a state-space dynamical system model for chemical thermodynamics. In particular, we use the law of mass-action to obtain the dynamics of chemical reaction networks. Furthermore, using the notion of the chemical potential [8,9], we unify our state space mass-action kinetics model with our thermodynamic dynamical system model involving energy exchange. In addition, we show that entropy production during chemical reactions is nonnegative and the dynamical system states of our chemical thermodynamic state space model converge to a state of temperature equipartition and zero affinity (i.e., the difference between the chemical potential of the reactants and the chemical potential of the products in a chemical reaction).

The central thesis of this paper is to present a state space formulation for equilibrium and nonequilibrium thermodynamics based on a dynamical system theory combined with interconnected nonlinear compartmental systems that ensures a consistent thermodynamic model for heat, energy, and mass flow. In particular, the proposed approach extends the framework developed in [4] addressing closed thermodynamic systems that exchange energy but not matter with the environment to open thermodynamic systems that exchange matter and energy with their environment. In addition, our results go beyond the results of [4] by developing rigorous notions of enthalpy, Gibbs free energy, Helmholtz free energy, and Gibbs’ chemical potential using a state space formulation of dynamics, energy and mass conservation principles, as well as the law of mass-action kinetics and the law of superposition of elementary reactions without invoking statistical mechanics arguments.

2. Notation, Definitions, and Mathematical Preliminaries

In this section, we establish notation, definitions, and provide some key results necessary for developing the main results of this paper. Specifically, $\mathbb{R}$ denotes the set of real numbers, ${\overline{\mathbb{Z}}}_{+}$ (respectively, ${\mathbb{Z}}_{+}$) denotes the set of nonnegative (respectively, positive) integers, ${\mathbb{R}}^q$ denotes the set of $q\times 1$ column vectors, ${\mathbb{R}}^{n\times m}$ denotes the set of $n\times m$ real matrices, ${\mathbb{P}}^n$ (respectively, ${\mathbb{N}}^n$) denotes the set of positive (respectively, nonnegative) definite matrices, ${(·)}^{\mathrm{T}}$ denotes transpose, $I_q$ or I denotes the $q\times q$ identity matrix, $\mathbf{e}$ denotes the ones vector of order q, that is, $\mathbf{e}\triangleq {[1,\dots ,1]}^{\mathrm{T}}\in {\mathbb{R}}^q$, and ${\mathbf{e}}_i\in {\mathbb{R}}^q$ denotes a vector with unity in the ith component and zeros elsewhere. For $x\in {\mathbb{R}}^q$ we write $x\ge \ge 0$ (respectively, $x>>0$) to indicate that every component of x is nonnegative (respectively, positive). In this case, we say that x is nonnegative or positive, respectively. Furthermore, ${\overline{\mathbb{R}}}_{+}^q$ and ${\mathbb{R}}_{+}^q$ denote the nonnegative and positive orthants of ${\mathbb{R}}^q$, that is, if $x\in {\mathbb{R}}^q$, then $x\in {\overline{\mathbb{R}}}_{+}^q$ and $x\in {\mathbb{R}}_{+}^q$ are equivalent, respectively, to $x\ge \ge 0$ and $x>>0$. Analogously, ${\overline{\mathbb{R}}}_{+}^{n\times m}$ (respectively, ${\mathbb{R}}_{+}^{n\times m}$) denotes the set of $n\times m$ real matrices whose entries are nonnegative (respectively, positive). For vectors $x,y\in {\mathbb{R}}^q$, with components $x_i$ and $y_i$, $i=1,\dots ,q$, we use $x\circ y$ to denote component-by-component multiplication, that is, $x\circ y\triangleq {[x_1{y}_1,\dots ,x_q{y}_q]}^{\mathrm{T}}$. Finally, we write $\partial \mathcal{S}$, $\stackrel{\circ }{\mathcal{S}}$, and $\overline{\mathcal{S}}$ to denote the boundary, the interior, and the closure of the set $\mathcal{S}$, respectively.

We write · for the Euclidean vector norm, $V^{\prime }\left(x\right)\triangleq \frac{\partial V\left(x\right)}{\partial x}$ for the Fréchet derivative of V at x, ${\mathcal{B}}_{\epsilon }\left(\alpha \right)$, $\alpha \in {\mathbb{R}}^q$, $\epsilon >0$, for the open ball centered at α with radius ε, and $x\left(t\right)\to \mathcal{M}$ as $t\to \infty$ to denote that $x\left(t\right)$ approaches the set $\mathcal{M}$ (that is, for every $\epsilon >0$ there exists $T>0$ such that dist$\left(x\right(t),\mathcal{M})<\epsilon$ for all $t>T$, where dist$(p,\mathcal{M})\triangleq {inf}_{x\in \mathcal{M}}\parallel p-x\parallel$). The notions of openness, convergence, continuity, and compactness that we use throughout the paper refer to the topology generated on $\mathcal{D}\subseteq {\mathbb{R}}^q$ by the norm $\parallel ·\parallel$. A subset $\mathcal{N}$ of $\mathcal{D}$ is relatively open in $\mathcal{D}$ if $\mathcal{N}$ is open in the subspace topology induced on $\mathcal{D}$ by the norm $\parallel ·\parallel$. A point $x\in {\mathbb{R}}^q$ is a subsequential limit of the sequence ${\left\{x_i\right\}}_{i=0}^{\infty }$ in ${\mathbb{R}}^q$ if there exists a subsequence of ${\left\{x_i\right\}}_{i=0}^{\infty }$ that converges to x in the norm $\parallel ·\parallel$. Recall that every bounded sequence has at least one subsequential limit. A divergent sequence is a sequence having no convergent subsequence.

Consider the nonlinear autonomous dynamical system

$$\begin{array}{ccc}\hfill \dot{x}\left(t\right)& =& f\left(x\left(t\right)\right),x\left(0\right)=x_0,t\in {\mathcal{I}}_{x_0}\hfill \end{array}$$

(1)

where $x\left(t\right)\in \subseteq {\mathbb{R}}^n$, $t\in {\mathcal{I}}_{x_0}$, is the system state vector, $\mathcal{D}$ is a relatively open set, $f:\mathcal{D}\to {\mathbb{R}}^n$ is continuous on $\mathcal{D}$, and ${\mathcal{I}}_{x_0}=[0,{\tau }_{x_0})$, $0\le {\tau }_{x_0}\le \infty$, is the maximal interval of existence for the solution $x(·)$ of Equation (1). We assume that, for every initial condition $x\left(0\right)\in \mathcal{D}$, the differential Equation (1) possesses a unique right-maximally defined continuously differentiable solution which is defined on $[0,\infty )$. Letting $s(·,x)$ denote the right-maximally defined solution of Equation (1) that satisfies the initial condition $x\left(0\right)=x$, the above assumptions imply that the map $s:[0,\infty )\times \mathcal{D}\to \mathcal{D}$ is continuous ([Theorem V.2.1] [10]), satisfies the consistency property $s(0,x)=x$, and possesses the semigroup property $s(t,s(\tau ,x\left)\right)=s(t+\tau ,x)$ for all $t,\tau \ge 0$ and $x\in \mathcal{D}$. Given $t\ge 0$ and $x\in \mathcal{D}$, we denote the map $s(t,·):\mathcal{D}\to \mathcal{D}$ by $s_t$ and the map $s(·,x):[0,\infty )\to \mathcal{D}$ by $s^x$. For every $t\in \mathbb{R}$, the map $s_t$ is a homeomorphism and has the inverse $s_{-t}$.

The orbit ${\mathcal{O}}_x$ of a point $x\in \mathcal{D}$ is the set $s^x\left([0,\infty )\right)$. A set ${\mathcal{D}}_{\mathrm{c}}\subseteq \mathcal{D}$ is positively invariant relative to Equation (1) if $s_t\left({\mathcal{D}}_{\mathrm{c}}\right)\subseteq {\mathcal{D}}_{\mathrm{c}}$ for all $t\ge 0$ or, equivalently, ${\mathcal{D}}_{\mathrm{c}}$ contains the orbits of all its points. The set ${\mathcal{D}}_{\mathrm{c}}$ is invariant relative to Equation (1) if $s_t\left({\mathcal{D}}_{\mathrm{c}}\right)={\mathcal{D}}_{\mathrm{c}}$ for all $t\ge 0$. The positive limit set of $x\in {\mathbb{R}}^q$ is the set $\omega \left(x\right)$ of all subsequential limits of sequences of the form ${\left\{s(t_i,x)\right\}}_{i=0}^{\infty }$, where ${\left\{t_i\right\}}_{i=0}^{\infty }$ is an increasing divergent sequence in $[0,\infty )$. $\omega \left(x\right)$ is closed and invariant, and ${\overline{\mathcal{O}}}_x={\mathcal{O}}_x\cup \omega \left(x\right)$ [7]. In addition, for every $x\in {\mathbb{R}}^q$ that has bounded positive orbits, $\omega \left(x\right)$ is nonempty and compact, and, for every neighborhood $\mathcal{N}$ of $\omega \left(x\right)$, there exists $T>0$ such that $s_t\left(x\right)\in \mathcal{N}$ for every $t>T$[7]. Furthermore, ∈ is an equilibrium point of Equation (1) if and only if $f\left(x_e\right)=0$ or, equivalently, $s(t,x_{\mathrm{e}})=x_{\mathrm{e}}$ for all $t\ge 0$. Finally, recall that if all solutions to Equation (1) are bounded, then it follows from the Peano–Cauchy theorem ([7] [p. 76]) that ${\mathcal{I}}_{x_0}=\mathbb{R}$.

Definition 2.1 ([11] [pp. 9, 10] ) Let $f={[f_1,\dots ,f_n]}^{\mathrm{T}}:\mathcal{D}\subseteq {\overline{\mathbb{R}}}_{+}^n\to {\mathbb{R}}^n$. Then f is essentially nonnegative if $f_i\left(x\right)\ge 0$, for all $i=1,\dots ,n$, and $x\in {\overline{\mathbb{R}}}_{+}^n$ such that $x_i=0$, where $x_i$ denotes the ith component of x.

Proposition 2.1 ([11] [p. 12] ) Suppose ${\overline{\mathbb{R}}}_{+}^n\subset \mathcal{D}$. Then ${\overline{\mathbb{R}}}_{+}^n$ is an invariant set with respect to Equation (1) if and only if $f:\mathcal{D}\to {\mathbb{R}}^n$ is essentially nonnegative.

Definition 2.2 ([11] [pp. 13, 23] ) An equilibrium solution $x\left(t\right)\equiv x_e\in {\overline{\mathbb{R}}}_{+}^n$ to Equation (1) is Lyapunov stable with respect to ${\overline{\mathbb{R}}}_{+}^n$ if, for all $\epsilon >0$, there exists $\delta =\delta \left(\epsilon \right)>0$ such that if $x\in {\mathcal{B}}_{\delta }\left(x_{\mathrm{e}}\right)\cap {\overline{\mathbb{R}}}_{+}^n$, then $x\left(t\right)\in {\mathcal{B}}_{\epsilon }\left(x_{\mathrm{e}}\right)\cap {\overline{\mathbb{R}}}_{+}^n$, $t\ge 0$. An equilibrium solution $x\left(t\right)\equiv x_e\in {\overline{\mathbb{R}}}_{+}^n$ to Equation (1) is semistable with respect to ${\overline{\mathbb{R}}}_{+}^n$ if it is Lyapunov stable with respect to ${\overline{\mathbb{R}}}_{+}^n$ and there exists δ > 0 such that if $x_0\in {\mathcal{B}}_{\delta }\left(x_e\right)\cap {\overline{\mathbb{R}}}_{+}^n$, then ${lim}_{t\to \infty }x\left(t\right)$ exists and corresponds to a Lyapunov stable equilibrium point with respect to ${\overline{\mathbb{R}}}_{+}^n$ The system given by Equation (1) is said to be semistable with respect to ${\overline{\mathbb{R}}}_{+}^n$if every equilibrium point of Equation (1) is semistable with respect to ${\overline{\mathbb{R}}}_{+}^n$ The system given by Equation (1) is said to be globally semistable with respect to ${\overline{\mathbb{R}}}_{+}^n$ if Equation (1) is semistable with respect to ${\overline{\mathbb{R}}}_{+}^n$ and, for every $x_0\in {\overline{\mathbb{R}}}_{+}^n,{lim}_{t\to \infty }x\left(t\right)$ exists.

Proposition 2.2 ([11] [p. 22]) Consider the nonlinear dynamical system given by Equation (1) where f is essentially nonnegative and let $x\in {\overline{\mathbb{R}}}_{+}^n$. If the positive limit set of Equation (1) contains a Lyapunov stable (with respect to ${\overline{\mathbb{R}}}_{+}^n$) equilibrium point y, then $y={lim}_{t\to \infty }s(t,x)$.

3. Interconnected Thermodynamic Systems: A State Space Energy Flow Perspective

The fundamental and unifying concept in the analysis of thermodynamic systems is the concept of energy. The energy of a state of a dynamical system is the measure of its ability to produce changes (motion) in its own system state as well as changes in the system states of its surroundings. These changes occur as a direct consequence of the energy flow between different subsystems within the dynamical system. Heat (energy) is a fundamental concept of thermodynamics involving the capacity of hot bodies (more energetic subsystems with higher energy gradients) to produce work. As in thermodynamic systems, dynamical systems can exhibit energy (due to friction) that becomes unavailable to do useful work. This in turn contributes to an increase in system entropy, a measure of the tendency of a system to lose the ability of performing useful work. In this section, we use the state space formalism to construct a mathematical model of a thermodynamic system that is consistent with basic thermodynamic principles.

Specifically, we consider a large-scale system model with a combination of subsystems (compartments or parts) that is perceived as a single entity. For each subsystem (compartment) making up the system, we postulate the existence of an energy state variable such that the knowledge of these subsystem state variables at any given time $t=t_0$, together with the knowledge of any inputs (heat fluxes) to each of the subsystems for time $t\ge t_0$, completely determines the behavior of the system for any given time $t\ge t_0$. Hence, the (energy) state of our dynamical system at time t is uniquely determined by the state at time $t_0$ and any external inputs for time $t\ge t_0$ and is independent of the state and inputs before time $t_0$.

More precisely, we consider a large-scale interconnected dynamical system composed of a large number of units with aggregated (or lumped) energy variables representing homogenous groups of these units. If all the units comprising the system are identical (that is, the system is perfectly homogeneous), then the behavior of the dynamical system can be captured by that of a single plenipotentiary unit. Alternatively, if every interacting system unit is distinct, then the resulting model constitutes a microscopic system. To develop a middle-ground thermodynamic model placed between complete aggregation (classical thermodynamics) and complete disaggregation (statistical thermodynamics), we subdivide the large-scale dynamical system into a finite number of compartments, each formed by a large number of homogeneous units. Each compartment represents the energy content of the different parts of the dynamical system, and different compartments interact by exchanging heat. Thus, our compartmental thermodynamic model utilizes subsystems or compartments with describe the energy distribution among distinct regions in space with intercompartmental flows representing the heat transfer between these regions. Decreasing the number of compartments results in a more aggregated or homogeneous model, whereas increasing the number of compartments leads to a higher degree of disaggregation resulting in a heterogeneous model.

To formulate our state space thermodynamic model, consider the interconnected dynamical system $\mathcal{G}$ shown in Figure 1 involving energy exchange between q interconnected subsystems. Let $E_i:[0,\infty )\to {\overline{\mathbb{R}}}_{+}$ denote the energy (and hence a nonnegative quantity) of the ith subsystem, let $S_i:[0,\infty )\to \mathbb{R}$ denote the external power (heat flux) supplied to (or extracted from) the ith subsystem, let ${\varphi }_{ij}:{\overline{\mathbb{R}}}_{+}^q\to \mathbb{R}$, $i\ne j,i,j=1,\dots ,q$, denote the net instantaneous rate of energy (heat) flow from the jth subsystem to the ith subsystem, and let ${\sigma }_{ii}:{\overline{\mathbb{R}}}_{+}^q\to {\overline{\mathbb{R}}}_{+},i=1,\dots ,q$, denote the instantaneous rate of energy (heat) dissipation from the ith subsystem to the environment. Here, we assume that ${\varphi }_{ij}:{\overline{\mathbb{R}}}_{+}^q\to \mathbb{R}$, $i\ne j$, $i,j=1,\dots ,q$, and ${\sigma }_{ii}:{\overline{\mathbb{R}}}_{+}^q\to {\overline{\mathbb{R}}}_{+}$, $i=1,\dots ,q$, are locally Lipschitz continuous on ${\overline{\mathbb{R}}}_{+}^q$ and $S_i:[0,\infty )\to \mathbb{R},i=1,\dots ,q$, are bounded piecewise continuous functions of time.

Figure 1. Interconnected dynamical system $\mathcal{G}$.

Figure 1. Interconnected dynamical system $\mathcal{G}$.

An energy balance for the ith subsystem yields

$$\begin{array}{ccc}\hfill E_i\left(T\right)& =& E_i\left(t_0\right)+\left[\sum _{j=1,j\ne i}^q{\int }_{t_0}^T{\varphi }_{ij}\left(E\left(t\right)\right)\mathrm{d}t\right]-{\int }_{t_0}^T{\sigma }_{ii}\left(E\left(t\right)\right)\mathrm{d}t+{\int }_{t_0}^T{S}_i\left(t\right)\mathrm{d}t,T\ge t_0\hfill \end{array}$$

(2)

or, equivalently, in vector form,

$$\begin{array}{ccc}\hfill E\left(T\right)& =& E\left(t_0\right)+{\int }_{t_0}^{T}w\left(E\left(t\right)\right)\mathrm{d}t-{\int }_{t_0}^{T}d\left(E\left(t\right)\right)\mathrm{d}t+{\int }_{t_0}^{T}S\left(t\right)\mathrm{d}t,T\ge t_0\hfill \end{array}$$

(3)

where $E\left(t\right)\triangleq {[E_1\left(t\right),\dots ,E_q\left(t\right)]}^{\mathrm{T}}$, $t\ge t_0$, is the system energy state, $d\left(E\left(t\right)\right)\triangleq [{\sigma }_{11}\left(E\left(t\right)\right),\dots ,$${\sigma }_{qq}{\left(E\left(t\right)\right)]}^{\mathrm{T}}$, $t\ge t_0$, is the system dissipation, $S\left(t\right)\triangleq {[S_1\left(t\right),\dots ,S_q\left(t\right)]}^{\mathrm{T}}$, $t\ge t_0$, is the system heat flux, and $w={[w_1,\dots ,w_q]}^{\mathrm{T}}:{\overline{\mathbb{R}}}_{+}^q\to {\mathbb{R}}^q$ is such that

$$\begin{array}{c}\hfill w_i\left(E\right)=\sum _{j=1,j\ne i}^q{\varphi }_{ij}\left(E\right),E\in {\overline{\mathbb{R}}}_{+}^q\end{array}$$

(4)

Since ${\varphi }_{ij}:{\overline{\mathbb{R}}}_{+}^q\to \mathbb{R}$, $i\ne j,i,j=1,\dots ,q$, denotes the net instantaneous rate of energy flow from the jth subsystem to the ith subsystem, it is clear that ${\varphi }_{ij}\left(E\right)=-{\varphi }_{ji}\left(E\right)$, $E\in {\overline{\mathbb{R}}}_{+}^q$, $i\ne j$, $i,j=1,\dots ,q$, which further implies that ${\mathbf{e}}^{\mathrm{T}}w\left(E\right)=0$, $E\in {\overline{\mathbb{R}}}_{+}^q$.

Note that Equation (2) yields a conservation of energy equation and implies that the energy stored in the ith subsystem is equal to the external energy supplied to (or extracted from) the ith subsystem plus the energy gained by the ith subsystem from all other subsystems due to subsystem coupling minus the energy dissipated from the ith subsystem to the environment. Equivalently, Equation (2) can be rewritten as

$$\begin{array}{c}\hfill {\dot{E}}_i\left(t\right)=\left[\sum _{j=1,j\ne i}^q{\varphi }_{ij}\left(E\left(t\right)\right)\right]-{\sigma }_{ii}\left(E\left(t\right)\right)+S_i\left(t\right),E_i\left(t_0\right)=E_{i0},t\ge t_0\end{array}$$

(5)

or, in vector form,

$$\begin{array}{ccc}\hfill \dot{E}\left(t\right)& =& w\left(E\left(t\right)\right)-d\left(E\left(t\right)\right)+S\left(t\right),E\left(t_0\right)=E_0,t\ge t_0\hfill \end{array}$$

(6)

where $E_0\triangleq {[E_{10},\dots ,E_{q0}]}^{\mathrm{T}}$, yielding a power balance equation that characterizes energy flow between subsystems of the interconnected dynamical system $\mathcal{G}$. We assume that ${\varphi }_{ij}\left(E\right)\ge 0,E\in {\overline{\mathbb{R}}}_{+}^q$, whenever $E_i=0$, $i\ne j$, $i,j=1,\dots ,q$, and ${\sigma }_{ii}\left(E\right)=0$, whenever $E_i=0$, $i=1,\dots ,q$. The above constraint implies that if the energy of the ith subsystem of $\mathcal{G}$ is zero, then this subsystem cannot supply any energy to its surroundings or dissipate energy to the environment. In this case, $w\left(E\right)-d\left(E\right),E\in {\overline{\mathbb{R}}}_{+}^q$, is essentially nonnegative [12]. Thus, if $S\left(t\right)\equiv 0$, then, by Proposition 2.1, the solutions to Equation (6) are nonnegative for all nonnegative initial conditions. See [4,11,12] for further details.

Since our thermodynamic compartmental model involves intercompartmental flows representing energy transfer between compartments, we can use graph-theoretic notions with undirected graph topologies (i.e., bidirectional energy flows) to capture the compartmental system interconnections. Graph theory [13,14] can be useful in the analysis of the connectivity properties of compartmental systems. In particular, an undirected graph can be constructed to capture a compartmental model in which the compartments are represented by nodes and the flows are represented by edges or arcs. In this case, the environment must also be considered as an additional node.

For the interconnected dynamical system $\mathcal{G}$ with the power balance Equation (6), we define a connectivity matrix $\mathcal{C}\in {\mathbb{R}}^{q\times q}$ such that for $i\ne j$, $i,j=1,\dots ,q$, ${\mathcal{C}}_{(i,j)}\triangleq 1$ if ${\varphi }_{ij}\left(E\right)\not\equiv 0$ and ${\mathcal{C}}_{(i,j)}\triangleq 0$ otherwise, and ${\mathcal{C}}_{(i,i)}\triangleq -{\sum }_{k=1,k\ne i}^q{\mathcal{C}}_{(k,i)}$, $i=1,\dots ,q$. (The negative of the connectivity matrix, that is, $-\mathcal{C}$, is known as the graph Laplacian in the literature.) Recall that if rank$\mathcal{C}=q-1$, then $\mathcal{G}$ is strongly connected [4] and energy exchange is possible between any two subsystems of $\mathcal{G}$.

The next definition introduces a notion of entropy for the interconnected dynamical system $\mathcal{G}$.

Definition 3.1 Consider the interconnected dynamical system $\mathcal{G}$ with the power balance Equation (6). A continuously differentiable, strictly concave function $\mathcal{S}:{\overline{\mathbb{R}}}_{+}^q\to \mathbb{R}$ is called the entropy function of $\mathcal{G}$ if

$$\begin{array}{c}\hfill \left(\frac{\partial \mathcal{S}\left(E\right)}{\partial E_i}-\frac{\partial \mathcal{S}\left(E\right)}{\partial E_j}\right){\varphi }_{ij}\left(E\right)\ge 0,E\in {\overline{\mathbb{R}}}_{+}^q,i\ne j,i,j=1,\dots ,q\end{array}$$

(7)

and $\frac{\partial \mathcal{S}\left(E\right)}{\partial E_i}=\frac{\partial \mathcal{S}\left(E\right)}{\partial E_j}$ if and only if ${\varphi }_{ij}\left(E\right)=0$ with ${\mathcal{C}}_{(i,j)}=1$, $i\ne j$, $i,j=1,\dots ,q$.

It follows from Definition 3.1 that for an isolated system $\mathcal{G}$, that is, $S\left(t\right)\equiv 0$ and $d\left(E\right)\equiv 0$, the entropy function of $\mathcal{G}$ is a nondecreasing function of time. To see this, note that

$$\begin{array}{ccc}\hfill \dot{\mathcal{S}}\left(E\right)& =& \frac{\partial \mathcal{S}\left(E\right)}{\partial E}\dot{E}\hfill \\ & =& \sum _{i=1}^q\frac{\partial \mathcal{S}\left(E\right)}{\partial E_i}\sum _{j=1,j\ne i}^q{\varphi }_{ij}\left(E\right)\hfill \\ & =& \sum _{i=1}^q\sum _{j=i+1}^q\left(\frac{\partial \mathcal{S}\left(E\right)}{\partial E_i}-\frac{\partial \mathcal{S}\left(E\right)}{\partial E_j}\right){\varphi }_{ij}\left(E\right)\hfill \\ & \ge & 0,E\in {\overline{\mathbb{R}}}_{+}^q\hfill \end{array}$$

(8)

where $\frac{\partial \mathcal{S}\left(E\right)}{\partial E}\triangleq \left[\frac{\partial \mathcal{S}\left(E\right)}{\partial E_1},\dots ,\frac{\partial \mathcal{S}\left(E\right)}{\partial E_q}\right]$ and where we used the fact that ${\varphi }_{ij}\left(E\right)=-{\varphi }_{ji}\left(E\right)$, $E\in {\overline{\mathbb{R}}}_{+}^q$, $i\ne j$, $i,j=1,\dots ,q$.

Proposition 3.1 Consider the isolated (i.e., $S\left(t\right)\equiv 0$ and $d\left(E\right)\equiv 0$) interconnected dynamical system $\mathcal{G}$ with the power balance Equation (6). Assume that rank$\mathcal{C}=q-1$ and there exists an entropy function $\mathcal{S}:{\overline{\mathbb{R}}}_{+}^q\to \mathbb{R}$ of $\mathcal{G}$. Then, ${\sum }_{j=1}^q{\varphi }_{ij}\left(E\right)=0$ for all $i=1,\dots ,q$ if and only if $\frac{\partial \mathcal{S}\left(E\right)}{\partial E_1}=\cdots =\frac{\partial \mathcal{S}\left(E\right)}{\partial E_q}$. Furthermore, the set of nonnegative equilibrium states of Equation (6) is given by ${\mathcal{E}}_0\triangleq \left\{E\in {\overline{\mathbb{R}}}_{+}^q:\frac{\partial \mathcal{S}\left(E\right)}{\partial E_1}=\cdots =\frac{\partial \mathcal{S}\left(E\right)}{\partial E_q}\right\}$.

Proof. If $\frac{\partial \mathcal{S}\left(E\right)}{\partial E_i}=\frac{\partial \mathcal{S}\left(E\right)}{\partial E_j}$, then ${\varphi }_{ij}\left(E\right)=0$ for all $i,j=1,\dots ,q$, which implies that ${\sum }_{j=1}^q{\varphi }_{ij}\left(E\right)=0$ for all $i=1,\dots ,q$. Conversely, assume that ${\sum }_{j=1}^q{\varphi }_{ij}\left(E\right)=0$ for all $i=1,\dots ,q$, and, since $\mathcal{S}$ is an entropy function of $\mathcal{G}$, it follows that

$$\begin{array}{ccc}\hfill 0& =& \sum _{i=1}^q\sum _{j=1}^q\frac{\partial \mathcal{S}\left(E\right)}{\partial E_i}{\varphi }_{ij}\left(E\right)\hfill \\ & =& \sum _{i=1}^{q-1}\sum _{j=i+1}^q\left(\frac{\partial \mathcal{S}\left(E\right)}{\partial E_i}-\frac{\partial \mathcal{S}\left(E\right)}{\partial E_j}\right){\varphi }_{ij}\left(E\right)\hfill \\ & \ge & 0\hfill \end{array}$$

where we have used the fact that ${\varphi }_{ij}\left(E\right)=-{\varphi }_{ji}\left(E\right)$ for all $i,j=1,\dots ,q$. Hence,

$$\left(\frac{\partial \mathcal{S}\left(E\right)}{\partial E_i}-\frac{\partial \mathcal{S}\left(E\right)}{\partial E_j}\right){\varphi }_{ij}\left(E\right)=0$$

for all $i,j=1,\dots ,q$. Now, the result follows from the fact that rank$\mathcal{C}=q-1$. □

Theorem 3.1 Consider the isolated (i.e., $S\left(t\right)\equiv 0$ and $d\left(E\right)\equiv 0$) interconnected dynamical system $\mathcal{G}$ with the power balance Equation (6). Assume that rank$\mathcal{C}=q-1$ and there exists an entropy function $\mathcal{S}:{\overline{\mathbb{R}}}_{+}^q\to \mathbb{R}$ of $\mathcal{G}$. Then the isolated system $\mathcal{G}$ is globally semistable with respect to ${\overline{\mathbb{R}}}_{+}^q$.

Proof. Since $w(·)$ is essentially nonnegative, it follows from Proposition 2.1 that $E\left(t\right)\in {\overline{\mathbb{R}}}_{+}^q$, $t\ge t_0$, for all $E_0\in {\overline{\mathbb{R}}}_{+}^q$. Furthermore, note that since ${\mathbf{e}}^{\mathrm{T}}w\left(E\right)=0$, $E\in {\overline{\mathbb{R}}}_{+}^q$, it follows that ${\mathbf{e}}^{\mathrm{T}}\dot{E}\left(t\right)=0$, $t\ge t_0$. In this case, ${\mathbf{e}}^{\mathrm{T}}E\left(t\right)={\mathbf{e}}^{\mathrm{T}}{E}_0$, $t\ge t_0$, which implies that $E\left(t\right)$, $t\ge t_0$, is bounded for all $E_0\in {\overline{\mathbb{R}}}_{+}^q$. Now, it follows from Equation (8) that $\mathcal{S}\left(E\right(t\left)\right)$, $t\ge t_0$, is a nondecreasing function of time, and hence, by the Krasovskii–LaSalle theorem [7], $E\left(t\right)\to \mathcal{R}\triangleq \{E\in {\overline{\mathbb{R}}}_{+}^q:\dot{\mathcal{S}}\left(E\right)=0\}$ as $t\to \infty$. Next, it follows from Equation (8), Definition 3.1, and the fact that rank$\mathcal{C}=q-1$, that $\mathcal{R}=\left\{E\in {\overline{\mathbb{R}}}_{+}^q:\frac{\partial \mathcal{S}\left(E\right)}{\partial E_1}=\cdots =\frac{\partial \mathcal{S}\left(E\right)}{\partial E_q}\right\}={\mathcal{E}}_0$.

Now, let $E_{\mathrm{e}}\in {\mathcal{E}}_0$ and consider the continuously differentiable function $V:{\mathbb{R}}^q\to \mathbb{R}$ defined by

$$V\left(E\right)\triangleq \mathcal{S}\left(E_{\mathrm{e}}\right)-\mathcal{S}\left(E\right)-{\lambda }_{\mathrm{e}}({\mathbf{e}}^{\mathrm{T}}{E}_{\mathrm{e}}-{\mathbf{e}}^{\mathrm{T}}E)$$

where ${\lambda }_{\mathrm{e}}\triangleq \frac{\partial \mathcal{S}}{\partial E_1}\left(E_{\mathrm{e}}\right)$. Next, note that $V\left(E_{\mathrm{e}}\right)=0$, $\frac{\partial V}{\partial E}\left(E_{\mathrm{e}}\right)=-\frac{\partial \mathcal{S}}{\partial E}\left(E_{\mathrm{e}}\right)+{\lambda }_{\mathrm{e}}{\mathbf{e}}^{\mathrm{T}}=0$, and, since $\mathcal{S}(·)$ is a strictly concave function, $\frac{{\partial }^{2}V}{\partial E^2}\left(E_{\mathrm{e}}\right)=-\frac{{\partial }^{2}S}{\partial E^2}\left(E_{\mathrm{e}}\right)>0$, which implies that $V(·)$ admits a local minimum at $E_{\mathrm{e}}$. Thus, $V\left(E_{\mathrm{e}}\right)=0$, there exists 0 such that $V\left(E\right)>0$, $E\in {\mathcal{B}}_{\delta }\left(E_{\mathrm{e}}\right)\backslash \left\{E_{\mathrm{e}}\right\}$, and $\dot{V}\left(E\right)=-\dot{\mathcal{S}}\left(E\right)\le 0$ for all $E\in {\mathcal{B}}_{\delta }\left(E_{\mathrm{e}}\right)\backslash \left\{E_{\mathrm{e}}\right\}$, which shows that $V(·)$ is a Lyapunov function for $\mathcal{G}$ and $E_{\mathrm{e}}$ is a Lyapunov stable equilibrium of $\mathcal{G}$. Finally, since, for every $E_0\in {\overline{\mathbb{R}}}_{+}^n$, $E\left(t\right)\to {\mathcal{E}}_0$ as $t\to \infty$ and every equilibrium point of $\mathcal{G}$ is Lyapunov stable, it follows from Proposition 2.2 that $\mathcal{G}$ is globally semistable with respect to ${\overline{\mathbb{R}}}_{+}^q$. □

In classical thermodynamics, the partial derivative of the system entropy with respect to the system energy defines the reciprocal of the system temperature. Thus, for the interconnected dynamical system $\mathcal{G}$,

$$\begin{array}{c}\hfill T_i\triangleq {\left(\frac{\partial \mathcal{S}\left(E\right)}{\partial E_i}\right)}^{-1},i=1,\dots ,q\end{array}$$

(9)

represents the temperature of the ith subsystem. Equation (7) is a manifestation of the second law of thermodynamics and implies that if the temperature of the jth subsystem is greater than the temperature of the ith subsystem, then energy (heat) flows from the jth subsystem to the ith subsystem. Furthermore, $\frac{\partial \mathcal{S}\left(E\right)}{\partial E_i}=\frac{\partial \mathcal{S}\left(E\right)}{\partial E_j}$ if and only if ${\varphi }_{ij}\left(E\right)=0$ with ${\mathcal{C}}_{(i,j)}=1$, $i\ne j$, $i,j=1,\dots ,q$, implies that temperature equality is a necessary and sufficient condition for thermal equilibrium. This is a statement of the zeroth law of thermodynamics. As a result, Theorem 3.1 shows that, for a strongly connected system $\mathcal{G}$, the subsystem energies converge to the set of equilibrium states where the temperatures of all subsystems are equal. This phenomenon is known as equipartition of temperature [4] and is an emergent behavior in thermodynamic systems. In particular, all the system energy is eventually transferred into heat at a uniform temperature, and hence, all dynamical processes in $\mathcal{G}$ (system motions) would cease.

The following result presents a sufficient condition for energy equipartition of the system, that is, the energies of all subsystems are equal. This state of energy equipartition is uniquely determined by the initial energy in the system.

Theorem 3.2 Consider the isolated (i.e., $S\left(t\right)\equiv 0$ and $d\left(E\right)\equiv 0$) interconnected dynamical system $\mathcal{G}$ with the power balance Equation (6). Assume that rank $\mathcal{C}=q-1$ and there exists a continuously differentiable, strictly concave function $f:{\overline{\mathbb{R}}}_{+}\to \mathbb{R}$ such that the entropy function $\mathcal{S}:{\overline{\mathbb{R}}}_{+}^q\to \mathbb{R}$ of $\mathcal{G}$ is given by $\mathcal{S}\left(E\right)={\sum }_{i=1}^{q}f\left(E_i\right)$. Then, the set of nonnegative equilibrium states of Equation (6) is given by ${\mathcal{E}}_0=\{\alpha \mathbf{e}:\alpha \ge 0\}$ and $\mathcal{G}$ is semistable with respect to ${\overline{\mathbb{R}}}_{+}^q$. Furthermore, $E\left(t\right)\to \frac{1}{q}\mathbf{e}{\mathbf{e}}^{\mathrm{T}}E\left(t_0\right)$ as $t\to \infty$ and $\frac{1}{q}\mathbf{e}{\mathbf{e}}^{\mathrm{T}}E\left(t_0\right)$ is a semistable equilibrium state of $\mathcal{G}$.

Proof. First, note that since $f(·)$ is a continuously differentiable, strictly concave function, it follows that

$$\left(\frac{\mathrm{d}f}{\mathrm{d}{E}_i}-\frac{\mathrm{d}f}{\mathrm{d}{E}_j}\right)(E_i-E_j)\le 0,E\in {\overline{\mathbb{R}}}_{+}^q,i,j=1,\dots ,q$$

which implies that Equation (7) is equivalent to

$$\left(E_i-E_j\right){\varphi }_{ij}\left(E\right)\le 0,E\in {\overline{\mathbb{R}}}_{+}^q,i\ne j,i,j=1,\dots ,q$$

and $E_i=E_j$ if and only if ${\varphi }_{ij}\left(E\right)=0$ with ${\mathcal{C}}_{(i,j)}=1$, $i\ne j$, $i,j=1,\dots ,q$. Hence, $-E^{\mathrm{T}}E$ is an entropy function of $\mathcal{G}$. Next, with $\mathcal{S}\left(E\right)=-\frac{1}{2}{E}^{\mathrm{T}}E$, it follows from Proposition 3.1 that ${\mathcal{E}}_0=\{\alpha \mathbf{e}\in {\overline{\mathbb{R}}}_{+}^q,\alpha \ge 0\}$. Now, it follows from Theorem 3.1 that $\mathcal{G}$ is globally semistable with respect to ${\overline{\mathbb{R}}}_{+}^q$. Finally, since ${\mathbf{e}}^{\mathrm{T}}E\left(t\right)={\mathbf{e}}^{\mathrm{T}}E\left(t_0\right)$ and $E\left(t\right)\to \mathcal{M}$ as $t\to \infty$, it follows that $E\left(t\right)\to \frac{1}{q}\mathbf{e}{\mathbf{e}}^{\mathrm{T}}E\left(t_0\right)$ as $t\to \infty$. Hence, with $\alpha =\frac{1}{q}{\mathbf{e}}^{\mathrm{T}}E\left(t_0\right)$, $\alpha \mathbf{e}=\frac{1}{q}\mathbf{e}{\mathbf{e}}^{\mathrm{T}}E\left(t_0\right)$ is a semistable equilibrium state of Equation (6). □

If $f\left(E_i\right)={log}_e(c+E_i)$, where $c>0$, so that $\mathcal{S}\left(E\right)={\sum }_{i=1}^q{log}_e(c+E_i)$, then it follows from Theorem 3.2 that ${\mathcal{E}}_0=\{\alpha \mathbf{e}:\alpha \ge 0\}$ and the isolated (i.e., $S\left(t\right)\equiv 0$ and $d\left(E\right)\equiv 0$) interconnected dynamical system $\mathcal{G}$ with the power balance Equation (6) is semistable. In this case, the absolute temperature of the ith compartment is given by $c+E_i$. Similarly, if $\mathcal{S}\left(E\right)=-\frac{1}{2}{E}^{\mathrm{T}}E$, then it follows from Theorem 3.2 that ${\mathcal{E}}_0=\{\alpha \mathbf{e}:\alpha \ge 0\}$ and the isolated (i.e., $S\left(t\right)\equiv 0$ and $d\left(E\right)\equiv 0$) interconnected dynamical system $\mathcal{G}$ with the power balance Equation (6) is semistable. In both cases, $E\left(t\right)\to \frac{1}{q}\mathbf{e}{\mathbf{e}}^{\mathrm{T}}E\left(t_0\right)$ as $t\to \infty$. This shows that the steady-state energy of the isolated interconnected dynamical system $\mathcal{G}$ is given by $\frac{1}{q}\mathbf{e}{\mathbf{e}}^{\mathrm{T}}E\left(t_0\right)=\frac{1}{q}{\sum }_{i=1}^q{E}_i\left(t_0\right)\mathbf{e}$, and hence is uniformly distributed over all subsystems of $\mathcal{G}$. This phenomenon is known as energy equipartition [4]. The aforementioned forms of $\mathcal{S}\left(E\right)$ were extensively discussed in the recent book [4] where $\mathcal{S}\left(E\right)={\sum }_{i=1}^q{log}_e(c+E_i)$ and $-\mathcal{S}\left(E\right)=\frac{1}{2}{E}^{\mathrm{T}}E$ are referred to, respectively, as the entropy and the ectropy functions of the interconnected dynamical system $\mathcal{G}$.

4. Work Energy, Gibbs Free Energy, Helmoholtz Free Energy, Enthalpy, and Entropy

In this section, we augment our thermodynamic energy flow model $\mathcal{G}$ with an additional (deformation) state representing subsystem volumes in order to introduce the notion of work into our thermodynamically consistent state space energy flow model. Specifically, we assume that each subsystem can perform (positive) work on the environment and the environment can perform (negative) work on the subsystems. The rate of work done by the ith subsystem on the environment is denoted by $d_{\mathrm{w}i}:{\overline{\mathbb{R}}}_{+}^q\times {\mathbb{R}}_{+}^q\to {\overline{\mathbb{R}}}_{+}$, $i=1,\dots ,q$, the rate of work done by the environment on the ith subsystem is denoted by $S_{\mathrm{w}i}:[0,\infty )\to {\overline{\mathbb{R}}}_{+}$, $i=1,\dots ,q$, and the volume of the ith subsystem is denoted by $V_i:[0,\infty )\to {\mathbb{R}}_{+}$, $i=1,\dots ,q$. The net work done by each subsystem on the environment satisfies

$$\begin{array}{c}\hfill p_i(E,V)\mathrm{d}{V}_i=(d_{\mathrm{w}i}(E,V)-S_{\mathrm{w}i}\left(t\right))\mathrm{d}t\end{array}$$

(10)

where $p_i(E,V)$, $i=1,\dots ,q$, denotes the pressure in the ith subsystem and $V\triangleq {[V_1,\dots ,V_q]}^{\mathrm{T}}$.

Furthermore, in the presence of work, the energy balance Equation (5) for each subsystem can be rewritten as

$$\begin{array}{c}\hfill \mathrm{d}{E}_i=w_i(E,V)\mathrm{d}t-(d_{\mathrm{w}i}(E,V)-S_{\mathrm{w}i}\left(t\right))\mathrm{d}t-{\sigma }_{ii}(E,V)\mathrm{d}t+S_i\left(t\right)\mathrm{d}t\end{array}$$

(11)

where $w_i(E,V)\triangleq {\sum }_{j=1,j\ne i}^q{\varphi }_{ij}(E,V)$, ${\varphi }_{ij}:{\overline{\mathbb{R}}}_{+}^q\times {\mathbb{R}}_{+}^q\to \mathbb{R}$, $i\ne j,i,j=1,\dots ,q$, denotes the net instantaneous rate of energy (heat) flow from the jth subsystem to the ith subsystem, ${\sigma }_{ii}:{\overline{\mathbb{R}}}_{+}^q\times {\mathbb{R}}_{+}^q\to {\overline{\mathbb{R}}}_{+}$, $i=1,\dots ,q$, denotes the instantaneous rate of energy dissipation from the ith subsystem to the environment, and, as in Section 3, $S_i:[0,\infty )\to \mathbb{R}$, $i=1,\dots ,q$, denotes the external power supplied to (or extracted from) the ith subsystem. It follows from Equations (10) and (11) that positive work done by a subsystem on the environment leads to a decrease in the internal energy of the subsystem and an increase in the subsystem volume, which is consistent with the first law of thermodynamics.

The definition of entropy for $\mathcal{G}$ in the presence of work remains the same as in Definition 3.1 with $\mathcal{S}\left(E\right)$ replaced by $\mathcal{S}(E,V)$ and with all other conditions in the definition holding for every $V>>0$. Next, consider the ith subsystem of $\mathcal{G}$ and assume that $E_j$ and $V_j$, $j\ne i$, $i=1,\dots ,q$, are constant. In this case, note that

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}\mathcal{S}}{\mathrm{d}t}& =& \frac{\partial \mathcal{S}}{\partial E_i}\frac{\mathrm{d}{E}_i}{\mathrm{d}t}+\frac{\partial \mathcal{S}}{\partial V_i}\frac{\mathrm{d}{V}_i}{\mathrm{d}t}\hfill \end{array}$$

(12)

and

$$\begin{array}{c}\hfill p_i(E,V)={\left(\frac{\partial \mathcal{S}}{\partial E_i}\right)}^{-1}\left(\frac{\partial \mathcal{S}}{\partial V_i}\right),i=1,\dots ,q\end{array}$$

(13)

It follows from Equations (10) and (11) that, in the presence of work energy, the power balance Equation (6) takes the new form involving energy and deformation states

$$\begin{array}{c}\hfill \dot{E}\left(t\right)=w(E\left(t\right),V\left(t\right))-d_{\mathrm{w}}(E\left(t\right),V\left(t\right))+S_{\mathrm{w}}\left(t\right)-d(E\left(t\right),V\left(t\right))+S\left(t\right),\\ \hfill E\left(t_0\right)=E_0,t\ge t_0,\end{array}$$

(14)

$$\begin{array}{ccc}\hfill \dot{V}\left(t\right)& =& D(E\left(t\right),V\left(t\right))(d_{\mathrm{w}}(E\left(t\right),V\left(t\right))-S_{\mathrm{w}}\left(t\right)),V\left(t_0\right)=V_0\hfill \end{array}$$

(15)

where $w(E,V)\triangleq {[w_1(E,V),\dots ,w_q(E,V)]}^{\mathrm{T}}$, $d_{\mathrm{w}}(E,V)\triangleq [d_{\mathrm{w}1}(E,V),\dots ,$ $d_{\mathrm{w}q}{(E,V)]}^{\mathrm{T}}$, $S_{\mathrm{w}}\left(t\right)$$\triangleq {[S_{\mathrm{w}1}\left(t\right),\dots ,S_{\mathrm{w}q}\left(t\right)]}^{\mathrm{T}}$, $d(E,V)\triangleq {[{\sigma }_{11}(E,V),\dots ,{\sigma }_{qq}(E,V)]}^{\mathrm{T}}$, $S\left(t\right)\triangleq {[S_1\left(t\right),\dots ,S_q\left(t\right)]}^{\mathrm{T}}$, and

$$\begin{array}{c}\hfill D(E,V)\triangleq \mathrm{diag}\left[\left(\frac{\partial \mathcal{S}}{\partial E_1}\right){\left(\frac{\partial \mathcal{S}}{\partial V_1}\right)}^{-1},\dots ,\left(\frac{\partial \mathcal{S}}{\partial E_q}\right){\left(\frac{\partial \mathcal{S}}{\partial V_q}\right)}^{-1}\right]\end{array}$$

(16)

Note that

$$\begin{array}{c}\hfill \left(\frac{\partial \mathcal{S}(E,V)}{\partial V}\right)D(E,V)=\frac{\partial \mathcal{S}(E,V)}{\partial E}\end{array}$$

(17)

The power balance and deformation Equations (14) and (15) represent a statement of the first law of thermodynamics. To see this, define the work L done by the interconnected dynamical system $\mathcal{G}$ over the time interval $[t_1,t_2]$ by

$$\begin{array}{c}\hfill L\triangleq {\int }_{t_1}^{t_2}{\mathbf{e}}^{\mathrm{T}}[d_{\mathrm{w}}(E\left(t\right),V\left(t\right))-S_{\mathrm{w}}\left(t\right)]\mathrm{d}t\end{array}$$

(18)

where ${[E^{\mathrm{T}}\left(t\right),V^{\mathrm{T}}\left(t\right)]}^{\mathrm{T}}$, $t\ge t_0$, is the solution to Equations (14) and (15). Now, premultiplying Equation (14) by ${\mathbf{e}}^{\mathrm{T}}$ and using the fact that ${\mathbf{e}}^{\mathrm{T}}w(E,V)=0$, it follows that

$$\begin{array}{c}\hfill \Delta U=-L+Q\end{array}$$

(19)

where $\Delta U=U\left(t_2\right)-U\left(t_1\right)\triangleq {\mathbf{e}}^{\mathrm{T}}E\left(t_2\right)-{\mathbf{e}}^{\mathrm{T}}E\left(t_1\right)$ denotes the variation in the total energy of the interconnected system $\mathcal{G}$ over the time interval $[t_1,t_2]$ and

$$\begin{array}{c}\hfill Q\triangleq {\int }_{t_1}^{t_2}{\mathbf{e}}^{\mathrm{T}}[S\left(t\right)-d(E\left(t\right),V\left(t\right))]\mathrm{d}t\end{array}$$

(20)

denotes the net energy received by $\mathcal{G}$ in forms other than work.

This is a statement of the first law of thermodynamics for the interconnected dynamical system $\mathcal{G}$ and gives a precise formulation of the equivalence between work and heat. This establishes that heat and mechanical work are two different aspects of energy. Finally, note that Equation (15) is consistent with the classical thermodynamic equation for the rate of work done by the system $\mathcal{G}$ on the environment. To see this, note that Equation (15) can be equivalently written as

$$\begin{array}{c}\hfill \mathrm{d}L={\mathbf{e}}^{\mathrm{T}}{D}^{-1}(E,V)\mathrm{d}V\end{array}$$

(21)

which, for a single subsystem with volume V and pressure p, has the classical form

$$\begin{array}{c}\hfill \mathrm{d}L=p\mathrm{d}V\end{array}$$

(22)

It follows from Definition 3.1 and Equations (14)–(17) that the time derivative of the entropy function satisfies

$$\begin{array}{cc}\hfill \dot{\mathcal{S}}(E,V)& =\frac{\partial \mathcal{S}(E,V)}{\partial E}\dot{E}+\frac{\partial \mathcal{S}(E,V)}{\partial V}\dot{V}\hfill \\ & =\frac{\partial \mathcal{S}(E,V)}{\partial E}w(E,V)-\frac{\partial \mathcal{S}(E,V)}{\partial E}(d_{\mathrm{w}}(E,V)-S_{\mathrm{w}}\left(t\right))\hfill \\ & -\frac{\partial \mathcal{S}(E,V)}{\partial E}(d(E,V)-S\left(t\right))+\frac{\partial \mathcal{S}(E,V)}{\partial V}D(E,V)(d_{\mathrm{w}}(E,V)-S_{\mathrm{w}}\left(t\right))\hfill \\ & =\sum _{i=1}^q\frac{\partial \mathcal{S}(E,V)}{\partial E_i}\sum _{j=1,j\ne i}^q{\varphi }_{ij}(E,V)+\sum _{i=1}^q\frac{\partial \mathcal{S}(E,V)}{\partial E_i}(S_i\left(t\right)-d_i(E,V))\hfill \\ & =\sum _{i=1}^q\sum _{j=i+1}^q\left(\frac{\partial \mathcal{S}(E,V)}{\partial E_i}-\frac{\partial \mathcal{S}(E,V)}{\partial E_j}\right){\varphi }_{ij}(E,V)\hfill \\ & +\sum _{i=1}^q\frac{\partial \mathcal{S}(E,V)}{\partial E_i}(S_i\left(t\right)-d_i(E,V))\hfill \\ & \ge \sum _{i=1}^q\frac{\partial \mathcal{S}(E,V)}{\partial E_i}(S_i\left(t\right)-d_i(E,V)),(E,V)\in {\overline{\mathbb{R}}}_{+}^q\times {\mathbb{R}}_{+}^q\hfill \end{array}$$

(23)

Noting that $\mathrm{d}{Q}_i\triangleq [S_i-{\sigma }_{ii}\left(E\right)]\mathrm{d}t$, $i=1,\dots ,q$, is the infinitesimal amount of the net heat received or dissipated by the ith subsystem of $\mathcal{G}$ over the infinitesimal time interval $\mathrm{d}t$, it follows from Equation (23) that

$$\begin{array}{c}\hfill \mathrm{d}\mathcal{S}\left(E\right)\ge \sum _{i=1}^q\frac{\mathrm{d}{Q}_i}{T_i}\end{array}$$

(24)

Inequality (24) is the classical Clausius inequality for the variation of entropy during an infinitesimal irreversible transformation.

Note that for an adiabatically isolated interconnected dynamical system (i.e., no heat exchange with the environment), Equation (23) yields the universal inequality

$$\begin{array}{c}\hfill \mathcal{S}(E\left(t_2\right),V\left(t_2\right))\ge \mathcal{S}(E\left(t_1\right),V\left(t_1\right)),t_2\ge t_1\end{array}$$

(25)

which implies that, for any dynamical change in an adiabatically isolated interconnected system $\mathcal{G}$, the entropy of the final system state can never be less than the entropy of the initial system state. In addition, in the case where $(E\left(t\right),V\left(t\right))\in {\mathcal{M}}_{\mathrm{e}}$, $t\ge t_0$, where ${\mathcal{M}}_{\mathrm{e}}\triangleq \{(E,V)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q:E=\alpha \mathbf{e},\alpha \ge 0,V\in {\mathbb{R}}_{+}^q\}$, it follows from Definition 3.1 and Equation (23) that Inequality (25) is satisfied as a strict inequality for all $(E,V)\in ({\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q)\backslash {\mathcal{M}}_{\mathrm{e}}$. Hence, it follows from Theorem 2.15 of [4] that the adiabatically isolated interconnected system $\mathcal{G}$ does not exhibit Poincaré recurrence in $({\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q)\backslash {\mathcal{M}}_{\mathrm{e}}$.

Next, we define the Gibbs free energy, the Helmholtz free energy, and the enthalpy functions for the interconnected dynamical system $\mathcal{G}$. For this exposition, we assume that the entropy of $\mathcal{G}$ is a sum of individual entropies of subsystems of $\mathcal{G}$, that is, $\mathcal{S}(E,V)={\sum }_{i=1}^q{\mathcal{S}}_i(E_i,V_i)$, $(E,V)\in {\overline{\mathbb{R}}}_{+}^q\times {\mathbb{R}}_{+}^q$. In this case, the Gibbs free energy of $\mathcal{G}$ is defined by

$$\begin{array}{ccc}\hfill G(E,V)& \triangleq & \hfill {\mathbf{e}}^{\mathrm{T}}E-\sum _{i=1}^q{\left(\frac{\partial \mathcal{S}(E,V)}{\partial E_i}\right)}^{-1}{\mathcal{S}}_i(E_i,V_i)+\sum _{i=1}^q{\left(\frac{\partial \mathcal{S}(E,V)}{\partial E_i}\right)}^{-1}\left(\frac{\partial \mathcal{S}(E,V)}{\partial V_i}\right)V_i\\ & & \hfill (E,V)\in {\overline{\mathbb{R}}}_{+}^q\times {\mathbb{R}}_{+}^q\end{array}$$

(26)

the Helmholtz free energy of $\mathcal{G}$ is defined by

$$\begin{array}{ccc}\hfill F(E,V)& \triangleq & {\mathbf{e}}^{\mathrm{T}}E-\sum _{i=1}^q{\left(\frac{\partial \mathcal{S}(E,V)}{\partial E_i}\right)}^{-1}{\mathcal{S}}_i(E_i,V_i),(E,V)\in {\overline{\mathbb{R}}}_{+}^q\times {\mathbb{R}}_{+}^q\hfill \end{array}$$

(27)

and the enthalpy of $\mathcal{G}$ is defined by

$$\begin{array}{ccc}\hfill H(E,V)& \triangleq & {\mathbf{e}}^{\mathrm{T}}E+\sum _{i=1}^q{\left(\frac{\partial \mathcal{S}(E,V)}{\partial E_i}\right)}^{-1}\left(\frac{\partial \mathcal{S}(E,V)}{\partial V_i}\right)V_i,(E,V)\in {\overline{\mathbb{R}}}_{+}^q\times {\mathbb{R}}_{+}^q\hfill \end{array}$$

(28)

Note that the above definitions for the Gibbs free energy, Helmholtz free energy, and enthalpy are consistent with the classical thermodynamic definitions given by $G(E,V)=U+pV-TS$, $F(E,V)=U-TS$, and $H(E,V)=U+pV$, respectively. Furthermore, note that if the interconnected system $\mathcal{G}$ is isothermal and isobaric, that is, the temperatures of subsystems of $\mathcal{G}$ are equal and remain constant with

$$\begin{array}{c}\hfill {\left(\frac{\partial \mathcal{S}(E,V)}{\partial E_1}\right)}^{-1}=\cdots ={\left(\frac{\partial \mathcal{S}(E,V)}{\partial E_q}\right)}^{-1}=T>0\end{array}$$

(29)

and the pressure $p_i(E,V)$ in each subsystem of $\mathcal{G}$ remains constant, respectively, then any transformation in $\mathcal{G}$ is reversible.

The time derivative of $G(E,V)$ along the trajectories of Equations (14) and (15) is given by

$$\begin{array}{ccc}\hfill \dot{G}(E,V)& =& {\mathbf{e}}^{\mathrm{T}}\dot{E}-\sum _{i=1}^q{\left(\frac{\partial \mathcal{S}(E,V)}{\partial E_i}\right)}^{-1}\left[\frac{\partial \mathcal{S}(E,V)}{\partial E_i}{\dot{E}}_i+\frac{\partial \mathcal{S}(E,V)}{\partial V_i}{\dot{V}}_i\right]\hfill \\ & & +\sum _{i=1}^q{\left(\frac{\partial \mathcal{S}(E,V)}{\partial E_i}\right)}^{-1}\left(\frac{\partial \mathcal{S}(E,V)}{\partial V_i}\right){\dot{V}}_i\hfill \\ & =& 0\hfill \end{array}$$

(30)

which is consistent with classical thermodynamics in the absence of chemical reactions.

For an isothermal interconnected dynamical system $\mathcal{G}$, the time derivative of $F(E,V)$ along the trajectories of Equations (14) and (15) is given by

$$\begin{array}{ccc}\hfill \dot{F}(E,V)& =& {\mathbf{e}}^{\mathrm{T}}\dot{E}-\sum _{i=1}^q{\left(\frac{\partial \mathcal{S}(E,V)}{\partial E_i}\right)}^{-1}\left[\frac{\partial \mathcal{S}(E,V)}{\partial E_i}{\dot{E}}_i+\frac{\partial \mathcal{S}(E,V)}{\partial V_i}{\dot{V}}_i\right]\hfill \\ & =& -\sum _{i=1}^q{\left(\frac{\partial \mathcal{S}(E,V)}{\partial E_i}\right)}^{-1}\left(\frac{\partial \mathcal{S}(E,V)}{\partial V_i}\right){\dot{V}}_i\hfill \\ & =& -\sum _{i=1}^q(d_{\mathrm{w}i}(E,V)-S_{\mathrm{w}i}\left(t\right))\hfill \\ & =& -L\hfill \end{array}$$

(31)

where L is the net amount of work done by the subsystems of $\mathcal{G}$ on the environment. Furthermore, note that if, in addition, the interconnected system $\mathcal{G}$ is isochoric, that is, the volumes of each of the subsystems of $\mathcal{G}$ remain constant, then $\dot{F}(E,V)=0$. As we see in the next section, in the presence of chemical reactions the interconnected system $\mathcal{G}$ evolves such that the Helmholtz free energy is minimized.

Finally, for the isolated ($S\left(t\right)\equiv 0$ and $d(E,V)\equiv 0$) interconnected dynamical system $\mathcal{G}$, the time derivative of $H(E,V)$ along the trajectories of Equations (14) and (15) is given by

$$\begin{array}{ccc}\hfill \dot{H}(E,V)& =& {\mathbf{e}}^{\mathrm{T}}\dot{E}+\sum _{i=1}^q{\left(\frac{\partial \mathcal{S}(E,V)}{\partial E_i}\right)}^{-1}\left(\frac{\partial \mathcal{S}(E,V)}{\partial V_i}\right){\dot{V}}_i\hfill \\ & =& {\mathbf{e}}^{\mathrm{T}}\dot{E}+\sum _{i=1}^q(d_{\mathrm{w}i}(E,V)-S_{\mathrm{w}i}\left(t\right))\hfill \\ & =& {\mathbf{e}}^{\mathrm{T}}w(E,V)\hfill \\ & =& 0\hfill \end{array}$$

(32)

5. Chemical Equilibria, Entropy Production, Chemical Potential, and Chemical Thermodynamics

In its most general form thermodynamics can also involve reacting mixtures and combustion. When a chemical reaction occurs, the bonds within molecules of the reactant are broken, and atoms and electrons rearrange to form products. The thermodynamic analysis of reactive systems can be addressed as an extension of the compartmental thermodynamic model described in Section 3 and Section 4. Specifically, in this case the compartments would qualitatively represent different quantities in the same space, and the intercompartmental flows would represent transformation rates in addition to transfer rates. In particular, the compartments would additionally represent quantities of different chemical substances contained within the compartment, and the compartmental flows would additionally characterize transformation rates of reactants into products. In this case, an additional mass balance is included for addressing conservation of energy as well as conservation of mass. This additional mass conservation equation would involve the law of mass-action enforcing proportionality between a particular reaction rate and the concentrations of the reactants, and the law of superposition of elementary reactions ensuring that the resultant rates for a particular species is the sum of the elementary reaction rates for the species.

In this section, we consider the interconnected dynamical system $\mathcal{G}$ where each subsystem represents a substance or species that can exchange energy with other substances as well as undergo chemical reactions with other substances forming products. Thus, the reactants and products of chemical reactions represent subsystems of $\mathcal{G}$ with the mechanisms of heat exchange between subsystems remaining the same as delineated in Section 3. Here, for simplicity of exposition, we do not consider work done by the subsystem on the environment or work done by the environment on the system. This extension can be easily addressed using the formulation in Section 4.

To develop a dynamical systems framework for thermodynamics with chemical reaction networks, let q be the total number of species (i.e., reactants and products), that is, the number of subsystems in $\mathcal{G}$, and let $X_j$, $j=1,\dots ,q$, denote the jth species. Consider a single chemical reaction described by

$$\begin{array}{c}\hfill \sum _{j=1}^q{A}_j{X}_j\stackrel{k}{⟶}\sum _{j=1}^q{B}_j{X}_j\end{array}$$

(33)

where $A_j$, $B_j$, $j=1,\dots ,q$, are the stoichiometric coefficients and k denotes the reaction rate. Note that the values of $A_j$ corresponding to the products and the values of $B_j$ corresponding to the reactants are zero. For example, for the familiar reaction

$$\begin{array}{c}\hfill 2{\mathrm{H}}_2+{\mathrm{O}}_2\stackrel{k}{⟶}2{\mathrm{H}}_2\mathrm{O}\end{array}$$

(34)

$X_1$, $X_2$, and $X_3$ denote the species ${\mathrm{H}}_2$, ${\mathrm{O}}_2$, and ${\mathrm{H}}_2\mathrm{O}$, respectively, and $A_1=2$, $A_2=1$, $A_3=0$, $B_1=0$, $B_2=0$, and $B_3=2$.

In general, for a reaction network consisting of $r\ge 1$ reactions, the ith reaction is written as

$$\begin{array}{c}\hfill \sum _{j=1}^q{A}_{ij}{X}_j\stackrel{k_i}{⟶}\sum _{j=1}^q{B}_{ij}{X}_j,i=1,\dots ,r\end{array}$$

(35)

where, for $i=1,\dots ,r$, $k_i>0$ is the reaction rate of the ith reaction, ${\sum }_{j=1}^q{A}_{ij}{X}_j$ is the reactant of the ith reaction, and ${\sum }_{j=1}^q{B}_{ij}{X}_j$ is the product of the ith reaction. Each stoichiometric coefficient $A_{ij}$ and $B_{ij}$ is a nonnegative integer. Note that each reaction in the reaction network given by Equation (35) is represented as being irreversible. Irreversibility here refers to the fact that part of the chemical reaction involves generation of products from the original reactants. Reversible chemical reactions that involve generation of products from the reactants and vice versa can be modeled as two irreversible reactions, one involving generation of products from the reactants and the other involving generation of the original reactants from the products. Hence, reversible reactions can be modeled by including the reverse reaction as a separate reaction. The reaction network given by Equation (35) can be written compactly in matrix-vector form as

$$\begin{array}{c}\hfill AX\stackrel{k}{⟶}BX\end{array}$$

(36)

where $X={[X_1,\dots ,X_q]}^{\mathrm{T}}$ is a column vector of species, $k={[k_1,\dots ,k_r]}^{\mathrm{T}}\in {\mathbb{R}}_{+}^r$ is a positive vector of reaction rates, and $A\in {\mathbb{R}}^{r\times q}$ and $B\in {\mathbb{R}}^{r\times q}$ are nonnegative matrices such that $A_{(i,j)}=A_{ij}$ and $B_{(i,j)}=B_{ij}$, $i=1,\dots ,r$, $j=1,\dots ,q$.

Let $n_j:[0,\infty )\to {\overline{\mathbb{R}}}_{+}$, $j=1,\dots ,q$, denote the mole number of the jth species and define $n\triangleq {[n_1,\dots ,n_q]}^{\mathrm{T}}$. Invoking the law of mass-action [15], which states that, for an elementary reaction, that is, a reaction in which all of the stoichiometric coefficients of the reactants are one, the rate of reaction is proportional to the product of the concentrations of the reactants, the species quantities change according to the dynamics [11,16]

$$\begin{array}{c}\hfill \dot{n}\left(t\right)={(B-A)}^{\mathrm{T}}K{n}^A\left(t\right),n\left(0\right)=n_0,t\ge t_0\end{array}$$

(37)

where $K\triangleq \mathrm{diag}[k_1,\dots ,k_r]\in {\mathbb{P}}^r$ and

$$\begin{array}{c}\hfill n^A\triangleq \left[\begin{array}{c}\prod _{j=1}^q{n}_j^{A_{1j}}\\ ⋮\\ \prod _{j=1}^q{n}_j^{A_{rj}}\end{array}\right]=\left[\begin{array}{c}{n}_1^{A_{11}}\cdots n_q^{A_{1q}}\\ ⋮\\ n_1^{A_{r1}}\cdots n_q^{A_{rq}}\end{array}\right]\in {\overline{\mathbb{R}}}_{+}^r\end{array}$$

(38)

For details regarding the law of mass-action and Equation (37), see [11,15,16,17]. Furthermore, let $M_j>0$, $j=1,\dots ,q$, denote the molar mass (i.e., the mass of one mole of a substance) of the jth species, let $m_j:[0,\infty )\to {\overline{\mathbb{R}}}_{+}$, $j=1,\dots ,q$, denote the mass of the jth species so that $m_j\left(t\right)=M_j{n}_j\left(t\right)$, $t\ge t_0$, $j=1,\dots ,q$, and let $m\triangleq {[m_1,\dots ,m_q]}^{\mathrm{T}}$. Then, using the transformation $m\left(t\right)=Mn\left(t\right)$, where $M\triangleq [M_1,\dots ,M_q]\in {\mathbb{P}}^q$, Equation (37) can be rewritten as the mass balance

$$\begin{array}{c}\hfill \dot{m}\left(t\right)=M{(B-A)}^{\mathrm{T}}\tilde{K}{m}^A\left(t\right),m\left(0\right)=m_0,t\ge t_0\end{array}$$

(39)

where $\tilde{K}\triangleq \left[\frac{k_1}{{\prod }_{j=1}^q{M}_j^{A_{1j}}},\dots ,\frac{k_r}{{\prod }_{j=1}^q{M}_j^{A_{rj}}}\right]\in {\mathbb{P}}^r$.

In the absence of nuclear reactions, the total mass of the species during each reaction in Equation (36) is conserved. Specifically, consider the ith reaction in Equation (36) given by Equation (35) where the mass of the reactants is ${\sum }_{j=1}^q{A}_{ij}{M}_j$ and the mass of the products is ${\sum }_{j=1}^q{B}_{ij}{M}_j$. Hence, conservation of mass in the ith reaction is characterized as

$$\begin{array}{c}\hfill \sum _{j=1}^q(B_{ij}-A_{ij})M_j=0,i=1,\dots ,r\end{array}$$

(40)

or, in general for Equation (36), as

$$\begin{array}{c}\hfill {\mathbf{e}}^{\mathrm{T}}M{(B-A)}^{\mathrm{T}}=0\end{array}$$

(41)

Note that it follows from Equations (39) and (41) that ${\mathbf{e}}^{\mathrm{T}}\dot{m}\left(t\right)\equiv 0$.

Equation (39) characterizes the change in masses of substances in the interconnected dynamical system $\mathcal{G}$ due to chemical reactions. In addition to the change of mass due to chemical reactions, each substance can exchange energy with other substances according to the energy flow mechanism described in Section 3; that is, energy flows from substances at a higher temperature to substances at a lower temperature. Furthermore, in the presence of chemical reactions, the exchange of matter affects the change of energy of each substance through the quantity known as the chemical potential.

The notion of the chemical potential was introduced by Gibbs in 1875–1878 [8,9] and goes far beyond the scope of chemistry, affecting virtually every process in nature [18,19,20]. The chemical potential has a strong connection with the second law of thermodynamics in that every process in nature evolves from a state of higher chemical potential towards a state of lower chemical potential. It was postulated by Gibbs [8,9] that the change in energy of a homogeneous substance is proportional to the change in mass of this substance with the coefficient of proportionality given by the chemical potential of the substance.

To elucidate this, assume the jth substance corresponds to the jth compartment and consider the rate of energy change of the jth substance of $\mathcal{G}$ in the presence of matter exchange. In this case, it follows from Equation (5) and Gibbs’ postulate that the rate of energy change of the jth substance is given by

$$\begin{array}{ccc}\hfill {\dot{E}}_j\left(t\right)& =& \hfill \left[\sum _{k=1,k\ne j}^q{\varphi }_{jk}\left(E\left(t\right)\right)\right]-{\sigma }_{jj}\left(E\left(t\right)\right)+S_j\left(t\right)+{\mu }_j(E\left(t\right),m\left(t\right)){\dot{m}}_j\left(t\right),E_j\left(t_0\right)=E_{j0},\\ & & \hfill t\ge t_0\end{array}$$

(42)

where ${\mu }_j:{\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q\to \mathbb{R}$, $j=1,\dots ,q$, is the chemical potential of the jth substance. It follows from Equation (42) that ${\mu }_j(·,·)$ is the chemical potential of a unit mass of the jth substance. We assume that if $E_j=0$, then ${\mu }_j(E,m)=0$, $j=1,\dots ,q$, which implies that if the energy of the jth substance is zero, then its chemical potential is also zero.

Next, using Equations (39) and (42), the energy and mass balances for the interconnected dynamical system $\mathcal{G}$ can be written as

$$\begin{array}{ccc}\hfill \dot{E}\left(t\right)& =& \hfill w\left(E\left(t\right)\right)+P(E\left(t\right),m\left(t\right))M{(B-A)}^{\mathrm{T}}\tilde{K}{m}^A\left(t\right)-d\left(E\left(t\right)\right)+S\left(t\right),E\left(t_0\right)=E_0,\\ & & \hfill t\ge t_0,\end{array}$$

(43)

$$\begin{array}{ccc}\hfill \dot{m}\left(t\right)& =& M{(B-A)}^{\mathrm{T}}\tilde{K}{m}^A\left(t\right),m\left(0\right)=m_0\hfill \end{array}$$

(44)

where $P(E,m)\triangleq [{\mu }_1(E,m),\dots ,{\mu }_q(E,m)]\in {\mathbb{R}}^{q\times q}$ and where $w(·)$, $d(·)$, and $S(·)$ are defined as in Section 3. It follows from Proposition 1 of [16] that the dynamics of Equation (44) are essentially nonnegative and, since ${\mu }_j(E,m)=0$ if $E_j=0$, $j=1,\dots ,q$, it also follows that, for the isolated dynamical system $\mathcal{G}$ (i.e., $S\left(t\right)\equiv 0$ and $d\left(E\right)\equiv 0$), the dynamics of Equations (43) and (44) are essentially nonnegative.

Note that, for the ith reaction in the reaction network given by Equation (36), the chemical potentials of the reactants and the products are ${\sum }_{j=1}^q{A}_{ij}{M}_j{\mu }_j(E,m)$ and ${\sum }_{j=1}^q{B}_{ij}{M}_j{\mu }_j(E,m)$, respectively. Thus,

$$\begin{array}{c}\hfill \sum _{j=1}^q{B}_{ij}{M}_j{\mu }_j(E,m)-\sum _{j=1}^q{A}_{ij}{M}_j{\mu }_j(E,m)\le 0,(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q\end{array}$$

(45)

is a restatement of the principle that a chemical reaction evolves from a state of a greater chemical potential to that of a lower chemical potential, which is consistent with the second law of thermodynamics. The difference between the chemical potential of the reactants and the chemical potential of the products is called affinity [21,22] and is given by

$$\begin{array}{c}\hfill {\nu }_i(E,m)=\sum _{j=1}^q{A}_{ij}{M}_j{\mu }_j(E,m)-\sum _{j=1}^q{B}_{ij}{M}_j{\mu }_j(E,m)\ge 0,i=1,\dots ,r\end{array}$$

(46)

Affinity is a driving force for chemical reactions and is equal to zero at the state of chemical equilibrium. A nonzero affinity implies that the system in not in equilibrium and that chemical reactions will continue to occur until the system reaches an equilibrium characterized by zero affinity. The next assumption provides a general form for the inequalities (45) and (46).

Assumption 5.1 For the chemical reaction network (36) with the mass balance Equation (44), assume that $\mu (E,m)>>0$ for all $E\ne 0$ and

$$\begin{array}{c}\hfill (B-A)M\mu (E,m)\le \le 0,(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q\end{array}$$

(47)

or, equivalently,

$$\begin{array}{c}\hfill \nu (E,m)=(A-B)M\mu (E,m)\ge \ge 0,(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q\end{array}$$

(48)

where $\mu (E,m)\triangleq {[{\mu }_1(E,m),\dots ,{\mu }_q(E,m)]}^{\mathrm{T}}$ is the vector of chemical potentials of the substances of $\mathcal{G}$ and $\nu (E,m)\triangleq {[{\nu }_1(E,m),\dots ,{\nu }_r(E,m)]}^{\mathrm{T}}$ is the affinity vector for the reaction network given by Equation (36).

Note that equality in Equation (47) or, equivalently, in Equation (48) characterizes the state of chemical equilibrium when the chemical potentials of the products and reactants are equal or, equivalently, when the affinity of each reaction is equal to zero. In this case, no reaction occurs and $\dot{m}\left(t\right)=0$, $t\ge t_0$.

Next, we characterize the entropy function for the interconnected dynamical system $\mathcal{G}$ with the energy and mass balances given by Equations (43) and (44). The definition of entropy for $\mathcal{G}$ in the presence of chemical reactions remains the same as in Definition 3.1 with $\mathcal{S}\left(E\right)$ replaced by $\mathcal{S}(E,m)$ and with all other conditions in the definition holding for every $m>>0$. Consider the jth subsystem of $\mathcal{G}$ and assume that $E_k$ and $m_k$, $k\ne j$, $k=1,\dots ,q$, are constant. In this case, note that

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}\mathcal{S}}{\mathrm{d}t}& =& \frac{\partial \mathcal{S}}{\partial E_j}\frac{\mathrm{d}{E}_j}{\mathrm{d}t}+\frac{\partial \mathcal{S}}{\partial m_j}\frac{\mathrm{d}{m}_j}{\mathrm{d}t}\hfill \end{array}$$

(49)

and recall that

$$\begin{array}{c}\hfill \frac{\partial \mathcal{S}}{\partial E}P(E,m)+\frac{\partial \mathcal{S}}{\partial m}=0\end{array}$$

(50)

Next, it follows from Equation (50) that the time derivative of the entropy function $\mathcal{S}(E,m)$ along the trajectories of Equations (43) and (44) is given by

$$\begin{array}{cc}\hfill \dot{\mathcal{S}}(E,m)& =\frac{\partial \mathcal{S}(E,m)}{\partial E}\dot{E}+\frac{\partial \mathcal{S}(E,m)}{\partial m}\dot{m}\hfill \\ & =\frac{\partial \mathcal{S}(E,m)}{\partial E}w\left(E\right)+\left(\frac{\partial \mathcal{S}(E,m)}{\partial E}P(E,m)+\frac{\partial \mathcal{S}(E,m)}{\partial m}\right)M{(B-A)}^{\mathrm{T}}\tilde{K}{m}^A\hfill \\ & +\frac{\partial \mathcal{S}(E,m)}{\partial E}S\left(t\right)-\frac{\partial \mathcal{S}(E,m)}{\partial E}d\left(E\right)\hfill \\ & =\frac{\partial \mathcal{S}(E,m)}{\partial E}w\left(E\right)+\frac{\partial \mathcal{S}(E,m)}{\partial E}S\left(t\right)-\frac{\partial \mathcal{S}(E,m)}{\partial E}d\left(E\right)\hfill \\ & =\sum _{i=1}^q\sum _{j=i+1}^q\left(\frac{\partial \mathcal{S}(E,m)}{\partial E_i}-\frac{\partial \mathcal{S}(E,m)}{\partial E_j}\right){\varphi }_{ij}\left(E\right)+\frac{\partial \mathcal{S}(E,m)}{\partial E}S\left(t\right)-\frac{\partial \mathcal{S}(E,m)}{\partial E}d\left(E\right),\hfill \end{array}$$

$$\begin{array}{cc}& (E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q\hfill \end{array}$$

(51)

For the isolated system $\mathcal{G}$ (i.e., $S\left(t\right)\equiv 0$ and $d\left(E\right)\equiv 0$), the entropy function of $\mathcal{G}$ is a nondecreasing function of time and, using identical arguments as in the proof of Theorem 3.1, it can be shown that $(E\left(t\right),m\left(t\right))\to \mathcal{R}\triangleq \left\{(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q:\frac{\partial \mathcal{S}(E,m)}{\partial E_1}=\cdots =\frac{\partial \mathcal{S}(E,m)}{\partial E_q}\right\}$ as $t\to \infty$ for all $(E_0,m_0)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q$.

The entropy production in the interconnected system $\mathcal{G}$ due to chemical reactions is given by

$$\begin{array}{ccc}\hfill \mathrm{d}{\mathcal{S}}_{\mathrm{i}}(E,m)& =& \frac{\partial \mathcal{S}(E,m)}{\partial m}\mathrm{d}m\hfill \\ & =& -\frac{\partial \mathcal{S}(E,m)}{\partial E}P(E,m)M{(B-A)}^{\mathrm{T}}\tilde{K}{m}^A\mathrm{d}t,(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q\hfill \end{array}$$

(52)

If the interconnected dynamical system $\mathcal{G}$ is isothermal, that is, all subsystems of $\mathcal{G}$ are at the same temperature

$$\begin{array}{c}\hfill {\left(\frac{\partial \mathcal{S}(E,m)}{\partial E_1}\right)}^{-1}=\cdots ={\left(\frac{\partial \mathcal{S}(E,m)}{\partial E_q}\right)}^{-1}=T\end{array}$$

(53)

where $T>0$ is the system temperature, then it follows from Assumption 5.1 that

$$\begin{array}{ccc}\hfill \mathrm{d}{\mathcal{S}}_{\mathrm{i}}(E,m)& =& -\frac{1}{T}{\mathbf{e}}^{\mathrm{T}}P(E,m)M{(B-A)}^{\mathrm{T}}\tilde{K}{m}^A\mathrm{d}t\hfill \\ & =& -\frac{1}{T}{\mu }^{\mathrm{T}}(E,m)M{(B-A)}^{\mathrm{T}}\tilde{K}{m}^A\mathrm{d}t\hfill \\ & =& \frac{1}{T}{\nu }^{\mathrm{T}}(E,m)\tilde{K}{m}^A\mathrm{d}t\hfill \\ & \ge & 0,(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q\hfill \end{array}$$

(54)

Note that since the affinity of a reaction is equal to zero at the state of a chemical equilibrium, it follows that equality in Equation (54) holds if and only if $\nu (E,m)=0$ for some $E\in {\overline{\mathbb{R}}}_{+}^q$ and $m\in {\overline{\mathbb{R}}}_{+}^q$.

Theorem 5.1 Consider the isolated (i.e., $S\left(t\right)\equiv 0$ and $d\left(E\right)\equiv 0$) interconnected dynamical system $\mathcal{G}$ with the power and mass balances given by Equations (43) and (44). Assume that rank $\mathcal{C}=q-1$, Assumption 5.1 holds, and there exists an entropy function $\mathcal{S}:{\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q\to \mathbb{R}$ of $\mathcal{G}$. Then $\left(E\right(t),m(t\left)\right)\to \mathcal{R}$ as $t\to \infty$, where $\left(E\right(t),m(t\left)\right)$, $t\ge t_0$, is the solution to Equations (43) and (44) with the initial condition $(E_0,m_0)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q$ and

$$\begin{array}{c}\hfill \mathcal{R}=\left\{(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q:\frac{\partial \mathcal{S}(E,m)}{\partial E_1}=\cdots =\frac{\partial \mathcal{S}(E,m)}{\partial E_q}\mathrm{and}\nu (E,m)=0\right\}\end{array}$$

(55)

where $\nu (·,·)$ is the affinity vector of $\mathcal{G}$.

Proof. Since the dynamics of the isolated system $\mathcal{G}$ are essentially nonnegative, it follows from Proposition 2.1 that $(E\left(t\right),m\left(t\right))\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q$, $t\ge t_0$, for all $(E_0,m_0)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q$. Consider a scalar function $v(E,m)={\mathbf{e}}^{\mathrm{T}}E+{\mathbf{e}}^{\mathrm{T}}m$, $(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q$, and note that $v(0,0)=0$ and $v(E,m)>0$, $(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q$, $(E,m)\ne (0,0)$. It follows from Equation (41), Assumption 5.1, and ${\mathbf{e}}^{\mathrm{T}}w\left(E\right)\equiv 0$ that the time derivative of $v(·,·)$ along the trajectories of Equations (43) and (44) satisfies

$$\begin{array}{ccc}\hfill \dot{v}(E,m)& =& {\mathbf{e}}^{\mathrm{T}}\dot{E}+{\mathbf{e}}^{\mathrm{T}}\dot{m}\hfill \\ & =& {\mathbf{e}}^{\mathrm{T}}P(E,m)M{(B-A)}^{\mathrm{T}}\tilde{K}{m}^A\hfill \\ & =& {\mu }^{\mathrm{T}}(E,m)M{(B-A)}^{\mathrm{T}}\tilde{K}{m}^A\hfill \\ & =& -{\nu }^{\mathrm{T}}(E,m)\tilde{K}{m}^A\hfill \\ & \le & 0,(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q\hfill \end{array}$$

(56)

which implies that the solution $\left(E\right(t),m(t\left)\right)$, $t\ge t_0$, to Equations (43) and (44) is bounded for all initial conditions $(E_0,m_0)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q$.

Next, consider the function $\tilde{v}(E,m)={\mathbf{e}}^{\mathrm{T}}E+{\mathbf{e}}^{\mathrm{T}}m-\mathcal{S}(E,m)$, $(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q$. Then it follows from Equations (51) and (56) that the time derivative of $\tilde{v}(·,·)$ along the trajectories of Equations (43) and (44) satisfies

$$\begin{array}{ccc}\hfill \dot{\tilde{v}}(E,m)& =& {\mathbf{e}}^{\mathrm{T}}\dot{E}+{\mathbf{e}}^{\mathrm{T}}\dot{m}-\dot{\mathcal{S}}(E,m)\hfill \\ & =& -{\nu }^{\mathrm{T}}(E,m)\tilde{K}{m}^A-\sum _{i=1}^q\sum _{j=i+1}^q\left(\frac{\partial \mathcal{S}(E,m)}{\partial E_i}-\frac{\partial \mathcal{S}(E,m)}{\partial E_j}\right){\varphi }_{ij}\left(E\right)\hfill \\ & \le & 0,(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q\hfill \end{array}$$

(57)

which implies that $\tilde{v}(·,·)$ is a nonincreasing function of time, and hence, by the Krasovskii–LaSalle theorem [7], $(E\left(t\right),m\left(t\right))\to \mathcal{R}\triangleq \{(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q:\dot{\tilde{v}}(E,m)=0\}$ as $t\to \infty$. Now, it follows from Definition 3.1, Assumption 5.1, and the fact that rank $\mathcal{C}=q-1$ that

$$\begin{array}{ccc}\hfill \mathcal{R}& =& \left\{(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q:\frac{\partial \mathcal{S}(E,m)}{\partial E_1}=\cdots =\frac{\partial \mathcal{S}(E,m)}{\partial E_q}\right\}\hfill \\ & & \cap \{(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q:\nu (E,m)=0\}\hfill \end{array}$$

(58)

which proves the result. □

Theorem 5.1 implies that the state of the interconnected dynamical system $\mathcal{G}$ converges to the state of thermal and chemical equilibrium when the temperatures of all substances of $\mathcal{G}$ are equal and the masses of all substances reach a state where all reaction affinities are zero corresponding to a halting of all chemical reactions.

Next, we assume that the entropy of the interconnected dynamical system $\mathcal{G}$ is a sum of individual entropies of subsystems of $\mathcal{G}$, that is, $\mathcal{S}(E,m)={\sum }_{j=1}^q{\mathcal{S}}_j(E_j,m_j)$, $(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q$. In this case, the Helmholtz free energy of $\mathcal{G}$ is given by

$$\begin{array}{c}\hfill F(E,m)={\mathbf{e}}^{\mathrm{T}}E-\sum _{j=1}^q{\left(\frac{\partial \mathcal{S}(E,m)}{\partial E_j}\right)}^{-1}{\mathcal{S}}_j(E_j,m_j),(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q\end{array}$$

(59)

If the interconnected dynamical system $\mathcal{G}$ is isothermal, then the derivative of $F(·,·)$ along the trajectories of Equations (43) and (44) is given by

$$\begin{array}{cc}\hfill \dot{F}(E,m)& ={\mathbf{e}}^{\mathrm{T}}\dot{E}-\sum _{j=1}^q{\left(\frac{\partial \mathcal{S}(E,m)}{\partial E_j}\right)}^{-1}{\dot{\mathcal{S}}}_j(E_j,m_j)\hfill \\ & ={\mathbf{e}}^{\mathrm{T}}\dot{E}-\sum _{j=1}^q{\left(\frac{\partial \mathcal{S}(E,m)}{\partial E_j}\right)}^{-1}\left[\frac{\partial {\mathcal{S}}_j(E_j,m_j)}{\partial E_j}{\dot{E}}_j+\frac{\partial {\mathcal{S}}_j(E_j,m_j)}{\partial m_j}{\dot{m}}_j\right]\hfill \\ & ={\mu }^{\mathrm{T}}(E,m)M{(B-A)}^{\mathrm{T}}\tilde{K}{m}^A\hfill \\ & =-{\nu }^{\mathrm{T}}(E,m)\tilde{K}{m}^A\hfill \\ & \le 0,(E,m)\in {\overline{\mathbb{R}}}_{+}^q\times {\overline{\mathbb{R}}}_{+}^q\hfill \end{array}$$

(60)

with equality in Equation (60) holding if and only if $\nu (E,m)=0$ for some $E\in {\overline{\mathbb{R}}}_{+}^q$ and $m\in {\overline{\mathbb{R}}}_{+}^q$, which determines the state of chemical equilibrium. Hence, the Helmholtz free energy of $\mathcal{G}$ evolves to a minimum when the pressure and temperature of each subsystem of $\mathcal{G}$ are maintained constant, which is consistent with classical thermodynamics. A similar conclusion can be arrived at for the Gibbs free energy if work energy considerations to and by the system are addressed. Thus, the Gibbs and Helmholtz free energies are a measure of the tendency for a reaction to take place in the interconnected system $\mathcal{G}$, and hence, provide a measure of the work done by the interconnected system $\mathcal{G}$.

6. Conclusion and Opportunities for Future Research

In this paper, we developed a system-theoretic perspective for classical thermodynamics and chemical reaction processes. In particular, we developed a nonlinear compartmental model involving heat flow, work energy, and chemical reactions that captures all of the key aspects of thermodynamics, including its fundamental laws. In addition, we showed that the interconnected compartmental model gives rise to globally semistable equilibria involving states of temperature equipartition. Finally, using the notion of the chemical potential, we combined our heat flow compartmental model with a state space mass-action kinetics model to capture energy and mass exchange in interconnected large-scale systems in the presence of chemical reactions. In this case, it was shown that the system states converge to a state of temperature equipartition and zero affinity.

The underlying intention of this paper as well as [4,5,6] has been to present one of the most useful and general physical branches of science in the language of dynamical systems theory. In particular, our goal has been to develop a dynamical system formalism of thermodynamics using a large-scale interconnected systems theory that bridges the gap between classical and statistical thermodynamics. The laws of thermodynamics are among the most firmly established laws of nature, and it is hoped that this work will help to stimulate increased interaction between physicists and dynamical systems and control theorists. Besides the fact that irreversible thermodynamics plays a critical role in the understanding of our physical universe, it forms the underpinning of several fundamental life science and engineering disciplines, including biological systems, physiological systems, neuroscience, chemical reaction systems, ecological systems, demographic systems, transportation systems, network systems, and power systems, to cite but a few examples.

An important area of science where the dynamical system framework of thermodynamics can prove invaluable is in neuroscience. Advances in neuroscience have been closely linked to mathematical modeling beginning with the integrate-and-fire model of Lapicque [23] and proceeding through the modeling of the action potential by Hodgkin and Huxley [24] to the current era of mathematical neuroscience; see [25,26] and the numerous references therein. Neuroscience has always had models to interpret experimental results from a high-level complex systems perspective; however, expressing these models with dynamic equations rather than words fosters precision, completeness, and self-consistency. Nonlinear dynamical system theory, and in particular system thermodynamics, is ideally suited for rigorously describing the behavior of large-scale networks of neurons.

Merging the two universalisms of thermodynamics and dynamical systems theory with neuroscience can provide the theoretical foundation for understanding the network properties of the brain by rigorously addressing large-scale interconnected biological neuronal network models that govern the neuroelectronic behavior of biological excitatory and inhibitory neuronal networks [27]. As in thermodynamics, neuroscience is a theory of large-scale systems wherein graph theory can be used in capturing the connectivity properties of system interconnections, with neurons represented by nodes, synapses represented by edges or arcs, and synaptic efficacy captured by edge weighting giving rise to a weighted adjacency matrix governing the underlying directed graph network topology. However, unlike thermodynamics, wherein energy spontaneously flows from a state of higher temperature to a state of lower temperature, neuron membrane potential variations occur due to ion species exchanges which evolve from regions of higher concentrations to regions of lower concentrations. And this evolution does not occur spontaneously but rather requires the opening and closing of specific gates within specific ion channels.

A particularly interesting application of nonlinear dynamical systems theory to the neurosciences is to study phenomena of the central nervous system that exhibit nearly discontinuous transitions between macroscopic states. A very challenging and clinically important problem exhibiting this phenomenon is the induction of general anesthesia [28,29,30,31,32]. In any specific patient, the transition from consciousness to unconsciousness as the concentration of anesthetic drugs increases is very sharp, resembling a thermodynamic phase transition. In current clinical practice of general anesthesia, potent drugs are administered which profoundly influence levels of consciousness and vital respiratory (ventilation and oxygenation) and cardiovascular (heart rate, blood pressure, and cardiac output) functions. These variation patterns of the physiologic parameters (i.e., ventilation, oxygenation, heart rate, blood pressure, and cardiac output) and their alteration with levels of consciousness can provide scale-invariant fractal temporal structures to characterize the degree of consciousness in sedated patients.

In particular, the degree of consciousness reflects the adaptability of the central nervous system and is proportional to the maximum work output under a fully conscious state divided by the work output of a given anesthetized state. A reduction in maximum work output (and oxygen consumption) or elevation in the anesthetized work output (or oxygen consumption) will thus reduce the degree of consciousness. Hence, the fractal nature (i.e., complexity) of conscious variability is a self-organizing emergent property of the large-scale interconnected biological neuronal network since it enables the central nervous system to maximize entropy production and optimally dissipate energy gradients. In physiologic terms, a fully conscious healthy patient would exhibit rich fractal patterns in space (e.g., fractal vasculature) and time (e.g., cardiopulmonary variability) that optimize the ability for oxygenation and ventilation. Within the context of aging and complexity in acute illnesses, variation of physiologic parameters and their relationship to system complexity, fractal variability, and system thermodynamics have been explored in [33,34,35,36,37,38].

Merging system thermodynamics with neuroscience can provide the theoretical foundation for understanding the mechanisms of action of general anesthesia using the network properties of the brain. Even though simplified mean field models have been extensively used in the mathematical neuroscience literature to describe large neural populations [26], complex large-scale interconnected systems are essential in identifying the mechanisms of action for general anesthesia [27]. Unconsciousness is associated with reduced physiologic parameter variability, which reflects the inability of the central nervous system to adopt, and thus, decomplexifying physiologic work cycles and decreasing energy consumption (ischemia, hypoxia) leading to a decrease in entropy production. The degree of consciousness is a function of the numerous coupling of the network properties in the brain that form a complex large-scale, interconnected system. Complexity here refers to the quality of a system wherein interacting subsystems self-organize to form hierarchical evolving structures exhibiting emergent system properties; hence, a complex dynamical system is a system that is greater than the sum of its subsystems or parts. This complex system—involving numerous nonlinear dynamical subsystem interactions making up the system—has inherent emergent properties that depend on the integrity of the entire dynamical system and not merely on a mean field simplified reduced-order model.

Developing a dynamical system framework for neuroscience [27] and merging it with system thermodynamics [4,5,6] by embedding thermodynamic state notions (i.e., entropy, energy, free energy, chemical potential, etc.) will allow us to directly address the otherwise mathematically complex and computationally prohibitive large-scale dynamical models that have been developed in the literature. In particular, a thermodynamically consistent neuroscience model would emulate the clinically observed self-organizing spatio-temporal fractal structures that optimally dissipate energy and optimize entropy production in thalamocortical circuits of fully conscious patients. This thermodynamically consistent neuroscience framework can provide the necessary tools involving semistability, synaptic drive equipartitioning (i.e., synchronization across time scales), energy dispersal, and entropy production for connecting biophysical findings to psychophysical phenomena for general anesthesia.

In particular, we conjecture that as the model dynamics transition to an aesthetic state the system will involve a reduction in system complexity—defined as a reduction in the degree of irregularity across time scales—exhibiting semistability and synchronization of neural oscillators (i.e., thermodynamic energy equipartitioning). In other words, unconsciousness will be characterized by system decomplexification. In addition, connections between thermodynamics, neuroscience, and the arrow of time [4,5,6] can be explored by developing an understanding of how the arrow of time is built into the very fabric of our conscious brain. Connections between thermodynamics and neuroscience is not limited to the study of consciousness in general anesthesia and can be seen in biochemical systems, ecosystems, gene regulation and cell replication, as well as numerous medical conditions (e.g., seizures, schizophrenia, hallucinations, etc.), which are obviously of great clinical importance but have been lacking rigorous theoretical frameworks. This is a subject of current research.