Port-Hamiltonian modeling of non-isothermal chemical reaction networks

Motivated by recent progress on the port-Hamiltonian formulation of isothermal chemical reaction networks and of the continuous stirred tank reactor, the present paper aims to develop a port-Hamiltonian formulation of chemical reaction networks in the non-isothermal case, and to exploit this for equilibrium and stability analysis.


Introduction
Modeling of chemical reaction networks has attracted much attention in the last decades due to its wide application in systems biology and chemical engineering.Previous work, such as [8,14,15], provides the foundation of a structural theory of isothermal chemical reaction networks governed by mass action kinetics.From then on, a series of papers about the modeling and analysis of mass action kinetics chemical reaction networks appeared (see [1,16,29,35]).In most of these papers, the chemical reaction networks are assumed to take place under isothermal condition.Consequently, the influence of in/outflow of heat can not be taken into account.Hence, non-isothermal chemical reaction networks still pose fundamental challenges.
In this paper, we aim to use the port-Hamiltonian framework for the modeling of non-isothermal mass action kinetics chemical reaction networks.Port-Hamiltonian systems theory (PHS) has been intensively employed in the modeling and passivitybased control of electrical, mechanical and electromechanical systems (see [20,30,31]).
In [32,33], a port-Hamiltonian formulation of isothermal mass action kinetics chemical reaction networks was provided.A first step to non-isothermal chemical reaction networks was taken in [36].Here, based on the previous works [5,7,24,25], a quasi port-Hamiltonian formulation for non-isothermal chemical reaction networks was developed.
The main contributions of the present paper are as follows.First, based on mass and energy balance equations, a port-Hamiltonian formulation for non-isothermal mass action kinetics chemical reaction networks which are detailed-balanced is developed.This formulation directly extends the port-Hamiltonian formulation of isothermal chemical reaction networks of [32,33], in contrast with the quasi port-Hamiltonian formulation in [36].It exhibits the energy balance and the thermodynamic principles in an explicit way.Based on the obtained port-Hamiltonian formulation, we provide a thermodynamic analysis of the existence and characterization of thermodynamic equilibria and their asymptotic stability.Being directly related with the energy and entropy functions, this port-Hamiltonian formulation is easily applicable to chemical and biological systems.The second contribution of this paper is the extension of the port-Hamiltonian formulation and the thermodynamic analysis to non-isothermal chemical reaction networks with external ports.
The structure of the paper is as follows.In Sect.2, some notation will be introduced which will be used in the remainder of the paper.Section 3 surveys the main elements of non-isothermal chemical reaction networks.Section 4 develops the port-Hamiltonian formulation of non-isothermal chemical reaction networks, and shows how this formulation is in line with the main laws of thermodynamics.In Sect.5, a thermodynamic analysis will be carried out, including the characterization of equilibria and their asymptotic stability.In Sect.6, an example-a genetic protein synthesis circuit with internal feedback and cell-to-cell communication-is discussed as an illustration of the developed theory.Section 7 extends the previous results to non-isothermal chemical reaction networks with external ports.

Notation
R m denotes the space of m-dimensional real vectors, and R m + the space of mdimensional real vectors whose entries are all strictly positive.The element-wise natural logarithm Ln : R m + → R m , x → Ln(x), is defined as the mapping whose ith component is given as (Ln(x)) i := ln(x i ).Similarly, Exp : R m + → R m , x → Exp(x), is the mapping whose ith component is given as (Exp(x)) i := exp(x i ).Note that Exp(x + z) = Exp(x)Exp(z), Ln(xz) = Ln(x) + Ln(z), and Ln( x z ) = Ln(x) − Ln(z), where xz ∈ R m is the element-wise product (xz) i := x i z i , i = 1, . . ., m, and x z ∈ R m is the element-wise quotient ( x z ) i = x i z i , i = 1, . . ., m.Also, we define the mapping Diag : R m → R m×m , v → Diag(v), where Diag(v) is the diagonal matrix with (Diag(v)) ii = v i .Finally, the i × i identity matrix is denoted as I i , the i × i zero matrix is denoted as 0 i×i , while the i × 1 zero vector is denoted as 0 i .The notation z tr is used to denote the transpose of the vector z.

The chemical reaction network structure
In this section, we will survey the basic topological structure of chemical reaction networks which will be used in the following sections.Consider a chemical reaction network composed of r reactions, m species and c complexes, given by the following reversible reaction scheme: with α i j , β i j being stoichiometric coefficients.The graph-theoretic formulation, according to [9,11,15], is to consider the chemical complexes defined by the left-hand and the right-hand sides of the chemical reactions, and to associate to each complex a vertex of a graph, while each reaction from left-hand to right-hand complex corresponds to a directed edge.The concentrations of the species are denoted as x i , i = 1, . . ., m, and the total vector of concentrations is denoted as x = [x 1 , . . ., x m ] tr .In order to capture the basic conservation laws of the chemical reactions, we define an m ×r matrix C, known as the stoichiometric matrix, whose (i, j)th element is the signed stoichiometric coefficient of the ith species in the jth reaction.Similarly, to define the connection between the complexes of each chemical reaction, we define an m × c matrix Z , called the complex stoichiometric matrix, whose ρth column captures the expression of the ρth complex in the ith chemical species.Any directed graph is characterized by an c × r matrix B, called the incidence matrix, whose (i, j)th element equals to −1 if vertex i is the tail vertex of edge j and 1 if vertex i is the head vertex of edge j , while 0 otherwise.The relation between these three matrices for the graph of complexes is The dynamics of the chemical reaction network can now be written as where v is the r-dimensional vector of reaction rates, whose jth element represents the jth reaction rate of the chemical reaction network.Each reaction is considered to be a combination of a forward reaction with forward rate equation and a reverse reaction with reverse rate equation, both given by mass action kinetics.Thus, the reaction rate of the jth reaction can be written as where Z i ρ is the (i, ρ)th element of the complex stoichiometric matrix Z .Here, k f j (T ), k b j (T ) are the forward/backward rate reaction coefficients of the jth reaction defined by the Arrhenius equation ( [4]) where E f j , E b j are the activation energies, k f j and k b j the forward and backward rate constants, and R is the Boltzmann constant.Using the element-wise natural logarithm Ln(x) defined in Sect.2, the jth element of rate equation vector can be equivalently written as where Z S j and Z P j denote the columns of the complex stoichiometry matrix Z corresponding to the substrate complex S j and the product complex P j of the jth reaction.

The standard form of balanced mass action and balanced energy action for non-isothermal chemical reaction networks
In this paper, we will focus on detailed balanced mass action kinetics chemical reaction networks (see [10,28,34]).First, we assume that in the chemical reaction network, the equilibration following any reaction event is much faster that any reaction time scale.Thus, all intensive thermodynamic variables are well defined and equal everywhere in the system.Then, we assume that the chemical reaction network is closed and undergoes an adiabatic process.That means there is no heat or mass transfer between the system and external environment.Moreover, the chemical reaction network is isochoric so that the volume change can be neglected, i.e. dV = 0.The definition of thermodynamic equilibrium for isothermal chemical reaction networks, see e.g.[32], is extended to non-isothermal networks as follows.Definition 4.1 A vector of concentrations x * is called an equilibrium for the dynamics ẋ = Cv(x, T ) for a certain temperature T if Cv(x * , T ) = 0, and a thermodynamic equilibrium if v(x * , T ) = 0.A chemical reaction network ẋ = Cv(x, T ) is called detailed-balanced if it admits a thermodynamic equilibrium for every temperature T .
In order to stress the dependence on T , the thermodynamical equilibrium will be denoted by x * (T ).The conditions for existence of a thermodynamic equilibrium will be discussed in Sect.5.1.Throughout this section we assume that there exists at least one thermodynamic equilibrium, like in the isothermal case, see e.g.[32,33].We will use the existence of this thermodynamic equilibrium to develop a port-Hamiltonian formulation.

Mass balance equations
Let us recall the mass balance equations of a detailed balanced reaction network according to [32].Let x * (T ) ∈ R m + be a thermodynamic equilibrium for a certain temperature T , i.e., v(x * (T ), T ) = 0 ( 7 ) Then we define the conductance κ j (T ) of the jth reaction as: Furthermore the reaction rate of the jth chemical reaction ( 6) can be rewritten as Now define the r × r diagonal matrix of conductances K (T ) as Collecting all the reaction rates in (9) and employing the incidence matrix B defined in Sect.3, the rate vector of a detailed balanced non-isothermal reaction network can be written as Hence the dynamics of a detailed balanced mass action kinetics reaction network can be expanded as where μ = RT Ln( x x * (T ) ) is the vector of chemical potentials and L := B K (T )B tr is the weighted Laplacian matrix for the reaction network graph, with weights given by the conductances κ 1 (T ), . . ., κ r (T ).
Note that the value of the conductances κ 1 (T ), . . ., κ r (T ) is not only dependent on the temperature T , but also on the choice of the thermodynamic equilibrium x * (T ).
However, if the reaction network graph is connected, then for any other thermodynamical equilibrium x * * (T ) for the same temperature T , there exists a positive constant d such that This property of the matrix K has been proved in [32].It implies that the dependence on x * (T ) is minor; choosing another thermodynamical equilibrium only involves a uniform scaling of K , and thus of L. Another well-known property of L is the fact that the matrix L is independent of the orientation of the graph [2].

Energy balance equations
In this section we will express the energy conservation for a closed chemical reaction network in order to encompass the thermodynamic properties of the system.
Assuming that in the system the variation of the volume may be neglected, i.e. dV = 0, Gibbs' relation reduces to where U denotes the internal energy, S the entropy, and the conjugated intensive variables are the chemical potential ∂U ∂ x = μ and the temperature ∂U ∂ S = T .This implies Using the Eq. ( 12), the first term on the right-hand side of ( 16) also equals Since the system is considered to be isolated, the energy balance equation is This implies that the second term in (17) equals In the next section we will combine these equations with (12) in order to derive a port-Hamiltonian formulation of non-isothermal and isolated reaction networks.123

Port-Hamiltonian formulation
In this section, we show how Sects.4.1 and 4.2 can be combined into a port-Hamiltonian formulation of the dynamics of detailed-balanced chemical reaction networks.Firstly, we define the state vector z = [ where x is the vector of concentrations and U the internal energy.Then we define the Hamiltonian function H = −S, where S is the entropy.Note that the Gibbs' relation (15) can also be written in the entropy formulation where d S dx i = − μ i T and d S dU = 1 T are the intensive thermodynamic variables conjugated to x i and the internal energy U .This implies that the co-state vector corresponding to Note that μ and T can be expressed as function of the components of this co-state vector.Now define the skew-symmetric matrix and the symmetric matrix It follows that the dynamics of the non-isothermal mass action kinetics chemical reaction network (1) given by the mass balance equation ( 12) and the energy balance equation (18), can be written into quasi port-Hamiltonian form As we will see in Sect.4.4, T μ tr Z LExp( Z tr μ RT ) ≥ 0 and thus R is positive semi-definite.The formulation ( 23) is called 'quasi port-Hamiltonian', since the structure matrices J and R depend on the co-state variables 123 port-Hamiltonian formulation.This formulation is comparable to the formulation of the mass balance and energy balance equations such as GENERIC, suggested in [18], or the port-Hamiltonian formulation with generating function being the availability function derived from the entropy function in [12].

Entropy balance equation
In this section, we shall relate the positive semi-definiteness of the dissipation matrix R in (22) with the second law of thermodynamics.With this in mind let us compute the time-derivative of the entropy RT .It has been shown in [28] (using the properties of the Laplacian matrix L) that for any γ ∈ R c while γ tr LExp(γ ) = 0 if and only if B tr γ = 0. Hence, the entropy balance equation becomes Here σ is the irreversible entropy source term.Note that in Eq. ( 26), the timederivative of the entropy S is deduced from the port-Hamiltonian formulation (23) defined in Sect.4.3.It is consistent with Eq. (19), which is deduced from the Gibbs' relation.
Furthermore, note that the positivity of the irreversible entropy source term is equivalent to the positive semi-definiteness of the dissipation matrix R in (22).Indeed the only non-zero term of R is the (m +1, m +1)th element, denoted as R m+1,m+1 , which is related to the entropy source term as In summary, the quasi port-Hamiltonian representation of chemical reaction networks given in ( 21), ( 22) and ( 23) represents the mass and energy balance equations.Moreover from its structure, it implies the entropy balance equation.It differs from the expression of energy and entropy balance equations in [23], which are expressed for the non-equilibrium biochemical systems and in [27], where the free energy and entropy balance are considered in an isothermal case when the species are diluted in a solvent, which acts as a thermal bath, while the pressure P is set by the environment.
It differs from the quasi port-Hamiltonian representation of the mass and entropy balance equations of chemical reaction networks in [26,36], by the fact that it is based on the energy balance equation instead of on the entropy balance equation.Note that the description based on the energy balance equation is classical [6], and more easily derived than the description based on the entropy balance equation.Moreover, the quasi-port-Hamiltonian formulation given in ( 21), ( 22) and ( 23) fundamentally differs from the representation of chemical reaction networks as port Hamiltonian systems in [21] as well, by the fact that this quasi port-Hamiltonian representation is established on the whole space of concentration vectors instead of only locally around an equilibrium point, as in [21].
Finally the quasi port-Hamiltonian directly extends the port-Hamiltonian formulation of isothermal chemical reaction networks obtained in [32,33] by including the energy balance equation.

Thermodynamic equilibria and asymptotic stability
The discussion in Sect. 4 is based on the assumption of existence of a thermodynamic equilibrium.Starting from the definition of a thermodynamic equilibrium of non-isothermal chemical reaction networks, we will derive in this section a full characterization of the set of equilibria, analogous to the case of isothermal chemical reaction networks in [32].
Subsequently, for stability analysis, we will use Lyapunov function as an availability function which is directly based on the quasi port-Hamiltonian representation given in (21), (22) and (23).Note that the use of availability functions for stability analysis is classical, see e.g.[12,13,19].

Thermodynamic equilibria
In this section, the existence of a thermodynamic equilibrium will be derived in the following linear-algebraic way [10].Recall the definition of a thermodynamic equilibrium for non-isothermal chemical reaction networks from Sect. 4. Let z * be a thermodynamic equilibrium under a certain temperature T , i.e., v(z * ) = 0.This implies that for any j = 1, . . ., r , 123 These equations are referred to as the detailed-balance equations.Denote K eq j (T ) Collecting all chemical reactions from 1 to r , and making use of the incidence matrix B, we obtain the following condition for a thermodynamical equilibrium x * (T ) where K eq is the r -dimensional vector with jth element K eq j , which is dependent on the temperature T .Therefore, for a given temperature T , there exists a thermodynamic equilibrium x * (T ) ∈ R m + if and only if k f j > 0, k b j > 0 for all j = 1, . . ., r , and In general, the equilibrium concentration x * (T ) may not be unique.Let x * * (T ) be another thermodynamic equilibrium for the same temperature T .Then That is to say, for a certain temperature T , once one thermodynamic equilibrium x * (T ) is given, the whole set of thermodynamic equilibria at the same temperature T can be found.Furthermore, since dU = 0, we have U * = U * * .Denote z * = (x * (T ), U * ) and z * * (T ) = (x * * (T ), U * * ), then it follows that the set of thermodynamic equilibria at the same temperature T can be written as This directly extends the classical result for isothermal chemical reaction networks, see e.g.[32].Note that the value of the terms Exp(C tr Ln(x * )) depend on temperature T , while the relation Exp(C tr Ln(x * * )) = Exp(C tr Ln(x * )) is not dependent on temperature T .
Since K eq (T ) = Exp(C tr Ln(x * )) as a function of T is monotone and injective, it follows that the set of thermodynamic equilibria

Asymptotic stability
For isothermal chemical reaction networks, it was shown in [15,32,33], that the Gibbs' free energy can be used as a Lyapunov function for proving asymptotic stability towards a unique equilibrium depending on the initial condition.In this section we aim at proving a similar result for the non-isothermal case based on the port-Hamiltonian formulation obtained in the previous section, employing the availability function.Note that this is different from [24,36], where an energy-based availability function was employed.
We define the entropy based availability function as where z o is a reference point taken as a thermodynamic equilibrium, cf.Sect. 4.
Theorem 5.1 Consider a detailed balanced chemical reaction network given by ( 21), ( 22) and (23), with A : R m+1 + → R given by (31).Then A has a strict minimum at z o with A(z o ) = 0, while the time-derivative d A dt is less than or equal to zero with equality only at z o .
Proof For homogeneous mixtures, the entropy function is necessarily concave [3].Moreover, the entropy is strict concave if at least one global extensive property (such as volume, total mass, or total mole number) is fixed [17].Recall the assumption that the chemical reaction network is isochoric, i.e. dV = 0, so the entropy is strict concave and A has a strict minimum at z o .Moreover, the time derivative of A(z) is given as where μ = RT Ln x x * is the vector of chemical potentials, x * .Since x * and x o are both thermodynamic equilibria, we obtain from Eq. ( 29) 123 Recall [28] that since L is a balanced weighted Laplacian matrix, for any γ ∈ R c , we have γ tr LExp(γ ) ≥ 0, while γ tr LExp(γ ) = 0 if and only if B tr γ = 0. Hence Therefore the time derivative of A(z) satisfies Let the system converge to a point denoted as and denote the equilibrium temperature associated with the equilibrium point z * by T * .We know that at equilibrium the entropy is maximal, implying that This is a classical statement in Chemical Engineering, and is comparable with the statement in [27], where for isothermal systems the Gibbs' free energy is minimized.According to Eq. ( 26), we have σ | z=z * = 0 (33) and for an isolated system we have By using the Eqs.( 33) and ( 34), the equilibrium point z * and T * can be determined.Then, by using a similar argument as in [11,32] , the following theorem will imply the asymptotic stability towards the set Σ T * .
Theorem 5.2 Consider the detailed-balanced chemical reaction network (21), ( 22) and ( 23) with T ∈ R + .Then for any + , where z * * (x * * , T * ) ∈ Σ T * is a thermodynamic equilibrium for temperature T * .As proved in [11,32], there exists a unique β ∈ ker C tr such that Combining with Theorem 4.1, this implies that the equilibrium z * is asymptotically stable with respect to all initial conditions in near z * .Hence the asymptotic stability of the quasi port-Hamiltonian system defined by ( 21), ( 22) and ( 23) is proved.

Example: a genetic circuit with internal feedback and cell-to-cell communication
The approach of the previous section will be illustrated on a chemical reaction network, taking place in a very common protein synthesis circuit in the cell of E. Coli in the large intestine of human beings [22].When the cell of E.Coli receives a 'message' from the environment (a kind of transcription process from extracellular space into the E. Coli cell), three chemical reactions will take place at the intercellular level: When the chemical reaction network reaches an equilibrium state, the cell will send out a 'message' to the environment (a reversed transcription process).This is a very efficient gene circuit for adjustment of the concentrations on different kind of protein in the E. Coli cell, with internal feedback and cell-to-cell communication.

Modeling
Let us denote the concentration of the species LuxR, AHL, LuxR − AHL, (LuxR − AHL) 2 , DNA and DNA − (LuxR − AHL) 2 as x 1 , . . ., x 6 .Hence, the state vector is defined as z = [x 1 , . . ., x 6 , U ] tr and the gradient vector of Hamiltonian function −S is given as d(−S) dz = [ μ 1 T , . . ., μ 6 T , 1 T ] tr .With m = 6, r = 3 and c = 5, the stoichiometric matrix C ∈ R 6×3 is written as Since the chemical reactions will take place naturally when the cell gets the 'message' from the environment, this means that the activation energies in Arrhenius equation are so small that they can be ignored, i.e., E f j = 0, j = 1, . . ., r and E b j = 0, j = 1, . . ., r .Hence, the matrix of conductances K becomes independent of T and takes the form which is independent of T .Therefore, the port-Hamiltonian formulation (23) for the genetic protein synthesis circuit is where the matrix J (−

Equilibrium and Lyapunov function
At thermodynamic equilibrium v(z * ) = 0, it can be verified that the equilibrium set Σ T is the 3-dimensional set given as To study its asymptotic stability, we define the availability function as in (31), where the reference point z o is taken to be a thermodynamic equilibrium under the temperature T .We have A(z) = 0 at z = z o , and as discussed in Sect.5.2, the time derivative of A(z) can be written as Therefore A(z) is a well-defined Lyapunov candidate.The port-Hamiltonian system (23) for the genetic protein synthesis circuit is asymptotically stable under temperature T .

Non-isothermal chemical reaction networks with ports
In many application areas, the chemical reaction networks under consideration are not isolated.That is to say, there exist mass exchange or heat exchange between the chemical reaction network and its environment.
In this section we will extend the port-Hamiltonian formulation for non-isothermal chemical reaction networks to the case of mass and heat exchange.As in the previous work (see [12,13,24,26]) when modeling and control of the Continuous Stirred Tank Reactor (CSTR), we define 'external ports' as the inflow/outflow of a mixture.Furthermore, we suppose that the output flow is such that the volume and pressure are constant [4].
Then Eq. ( 2) can be rewritten as where the vectors F e and F s are respectively the input and output concentration flows.The previous formulation ( 23) can be extended to the non-isothermal chemical reaction network with external ports as where as before the Hamiltonian function is given as H = −S.The internal energy can be written as where c pi , u 0i T 0 are respectively the heat capacity at constant pressure, reference molar energy and reference temperature.With constant volume and pressure the balance equation for the internal energy U can be written as where Q is the heat flux from the environment, and F ei and F si are the ith element of F e and F s .Furthermore, h ei and h si are respectively the input and output specific enthalpies.
Note that in the port-Hamiltonian formulation for non-isothermal chemical reaction network with ports (38), we still use the thermodynamic equilibrium z * (T ) for the chemical reaction network without ports under given temperature T , as defined in Sect. 4.
As before, in order to verify the asymptotic stability, we define the availability Eq. (31) as: It then remains to prove that A(z) is a Lyapunov function.Obviously, we have A(z o ) = 0.
Here, we assume that ΔF = F e − F s is a vector which can be described by mass action kinetics.That is to say, it can be developed as tr a constant balanced weighted Laplacian matrix corresponding to some incidence matrix B and ΔU = 0. Note that under this assumption, the external ports added to the network can be considered as another chemical reaction network, with the same complex Z , but with the different incidence matrix B and different rate coefficients K .Then the time derivative of A(z) becomes dt ≤ 0, and thus the port-Hamiltonian system for non-isothermal chemical reaction networks with ports is asymptotically stable for temperature T .Let us illustrate this on the following example.

Example: a genetic circuit with internal feedback and cell-to-cell communication (continued)
Recall the genetic protein synthesis circuit described in Sect. 1. Assume that there exists a port of flows ΔF = −Z L Exp(Z tr Ln x x * ) with   Hence A(z) is a well-defined Lyapunov candidate.We conclude that the port-Hamiltonian system for genetic protein synthesis circuit with a specific port ΔF is asymptotically stable for temperature T .

Conclusions and outlook
In this paper, a (quasi) port-Hamiltonian formulation has been developed for nonisothermal mass action kinetics chemical reaction networks.As an extensive result of the port-Hamiltonian formulation for isothermal chemical reaction network, and based on the mass balance and energy balance equations, this port-Hamiltonian formulation provides us a very explicit way to represent the chemical reaction networks and their thermodynamic properties, including the entropy balance and the conditions for the existence of thermodynamic equilibrium.As for the asymptotic stability, a comparable statement with the one in [27] has been found.
Moreover, this (quasi) port-Hamiltonian formulation and its thermodynamic analysis have been extended to non-isothermal chemical reaction networks with external ports.The results have been illustrated on a chemical reaction network in our body: the genetic circuit with internal feedback and cell-to-cell communication.
The focus of future work will be on the extension of current results to the modeling of interconnection of non-isothermal chemical reaction networks.Inspired by Rao and Esposito [27] and Qian and Beard [23], the interconnected port-Hamiltonian formulation will be considered as the combination of two driven (or chemostatted) chemical reaction networks with shared species.

123 where l 1 ,
l 2 and l 3 are positive constants.The port-Hamiltonian formulation (