An Observer for Mass-action Chemical Reaction Networks

European Journal of Control (2009)5:578–593 # 2009 EUCA DOI:10.3166/EJC.15.578–593

An Observer for Mass-action Chemical Reaction Networks Marcello Farina1,, Sergio Bittanti1, 1

Dipartimento di Elettronica e Informazione, Politecnico di Milano, Italy

Several are the challenges posed by the biological and biomedical research, especially for diagnostics and therapy of diseases. Among them a main task is the understanding of key cellular phenomena via the

comprehension of the interactions taking place in living organisms to form a unique complex. Initially the focus was on the taxonomic description of single species and subsystems. In the last decades the research has shifted towards the analysis of the dynamic interrelations among subsystems. This has been also spurred by the progresses in the fields of genomics, proteomics and measurement technologies, leading to availability of large amounts of data. Such developments paved the way to a novel discipline, systems biology. Research in systems biology takes advantage of a multidisciplinary approach, with emphasis on mathematical modelling and computer simulation [30, 35, 38, 39]. In this framework, understanding the organisms functioning on a cellular and subcellular scale should be achieved at different levels, [29], e.g. by unraveling the system’s structure, analyzing its dynamical behaviour and studying the intrinsic selfregulation mechanisms [37] of the cell. At the workplace of system dynamics study, experimental data analysis plays an important role since it enables quantitative understanding. Monitoring the level of biochemical species involved in interconnected reactions (e.g., metabolic reaction networks) is a major issue for the understanding of biochemical processes [1, 2, 32]. Once a reliable mathematical model of the considered biochemical process is worked out (see for example [5, 31]), the problem can be treated by resorting to the methodologies of systems and control theory [21, 36]. In particular, the problem of estimating the species concentration in the reaction networks from the

Correspondence to: S. Bittanti, E-mail: [email protected] E-mail: [email protected]

Received 19 May 2008; Accepted 4 May 2009 Recommended by D. Normand-Cyrot, S. Monaco

In biological research, experimental data analysis plays an important role since it enables quantitative understanding of biochemical processes. On the other hand, today’s measurement techniques, in continuous development, generally allow measuring a subset of the major system’s variables. Such major issue can be tackled by relying on a system’s and mathematical approach. For instance, first principles modelling of metabolic or signal transduction networks typically leads to a set of nonlinear differential equations. In this paper, we devise a nonlinear observer specifically suited for models of biochemical reaction networks. We show that the observer is locally convergent under certain observability conditions which can be inferred by elementary network analysis. The applicability and performance of the outlined observer are shown considering the state estimation problem for a benchmark biochemical reaction network. Keywords: Nonlinear observer design, Observability, Mass-action reaction networks, Systems biology

1. Introduction

579

An Observer for Mass-action Reaction Networks

available measurements can be tackled as the problem of designing a suitable (nonlinear) observer [27]. The inspiring paper [4] deals with the problem of state estimation for a subclass of the mass-action systems [13–15, 17, 22]. More general classes of massaction networks can be encountered, for example to encompass the effect of the crosstalk phenomena with different networks, or to describe formation of biochemical species. In this paper we consider a general class of weakly reversible mass-action networks (as in [9, 10]). The main contribution is the derivation of a suitable Luenberger state estimator, whose equations are inspired by the network structure. We will also show that the observer is locally convergent under certain observability conditions which can be inferred by elementary network analysis. In Section 2 we recall the major notions used to describe the topology and the dynamics of a general type mass-action reaction network. We also review some of the results in analysis of a mass-action network [14, 23–25, 34]. In Section 3 the nonlinear observer for the estimation of concentrations is described. To be precise, in Section 3.1 we work out the observer for a general mass-action network. Then, in Section 3.2 we apply the proposed observer to the subclass of the so-called zero deficiency networks. A result on the application of the observer in order to encompass degradation and formation reactions is provided in Section 3.3. Section 4 is devoted to a case study. Specifically, to the so-called Edelstein network investigated in [11] and originally presented in 7. Finally some concluding remarks are given in Section 5. For the ease of paper readability, all proofs are moved to the Appendix.

2. Models of Chemical Reaction Networks In the last decades we witnessed a growing interest towards models of mass-action chemical reaction networks [16, 26], also referred to as Horn–Jackson– Feinberg models [22]. In the 1970s the steady-state properties of such models have been successfully investigated [14, 23–25]. Furthermore, some major steady-state properties like complex balancing [12] have been introduced and studied. For an excellent discussion on such results see also [11]. Recently new research directions have been investigated [4, 34]. In order to provide a reliable background to the reader and introduce the basic definitions we first outline the basic concepts and main results available. Then in Section 2.2 we present a formalization of the input/ output dynamical model previously introduced by the

authors in [9, 10]. In Section 2.3 we recall the so-called zero deficiency theorem. Its main statement guarantees existence and uniqueness of steady state in a positive invariant space for the so-called zero deficiency networks. Finally, in Section 2.4 we account for formation and degradation reactions in a general mass-action network and we make some remarks on the analysis of the outcoming models. Some more advanced notions (which are needed in the technical proofs of the main statements of Section 3) are moved to the Appendix.

2.1. Topology of a Network The fundamental units considered in biology are the species. Species are chemicals, e.g. proteins, genes, sugars, small molecules or large compounds which can be detected at cellular and subcellular level. By complex we mean a set of species. A single reaction scheme comprises a reactant complex XR (composed by certain species) and a product complex XP (in general composed by other species). The species contained in the reactant complex are degraded by the reaction, and the species of the product complex are created. The stoichiometric coefficients indicate the proportion at which each species degrades (or is created) with respect to the whole complex degradation (or formation) rate. A reaction network is constituted by a set of reactions. We denote by n the total number of reacting chemicals and by  the number of complexes. To each complex we associate a vector bj ; j ¼ 1;    ; . The i-th element of bj is the stoichiometric amount of species i in the complex j. Thus, each vector bj is a n-dimensional vector (of integer entries), bj 2 Rn . The set of all vectors bj can be gathered in a unique matrix   B ¼ b 1 b2    b  ;

bj 2 R n ;

j ¼ 1;    ; 

ð1Þ

denominated complex matrix of dimension n  . Through this paper, we make the following assumption: Assumption 1: The matrix B has no zero rows. That is, the species considered in the reaction network model must be part of at least one complex, i.e. concentrations of species which neither react nor are the product of a given reaction should not be considered as state variables. We can describe each reaction in terms of the complexes involved (i.e., reactant and product complexes). For instance, consider the k-th reaction, transforming reactant (jR -th complex) into the

580

M. Farina and S. Bittanti

product (jP -th complex). Such reaction is characterized by the corresponding input and output stoichiometric vectors:

We define the matrices S; Sout ; Sin 2 Rr as: S ¼ ½ s1 s2    sr 

ð7aÞ

nink ¼ bjR

ð2aÞ

  r Sout ¼ s1out s2out    sout

ð7bÞ

k ¼ bjP nout

ð2bÞ

  Sin ¼ s1in s2in    srin

ð7cÞ

The net amount of each species involved in the k-th reaction is therefore nk ¼



k nout

ð2cÞ

nink

Letting r be the total number of reactions taking place in the network, we define the stoichiometric matrix, the output stoichiometric matrix and the input stoichiometric matrix (see [10]) as: N ¼ ½n1 n2    nr  ; nk 2 Rn ; j ¼ 1; :::; r   k 2 Rn ; Nout ¼ n1out n2out    nrout ; nout j ¼ 1; :::; r   Nin ¼ n1in n2in    nrin ; nink 2 Rn ; j ¼ 1; :::; r

ð4Þ

ð5Þ

sink ¼ wjR

ð6aÞ

k ¼ wjP sout

ð6bÞ

sk ¼



sink

1 Ð C2 2

3 C2 þ 2 C3 ! C4

Given a network comprising  different complexes, R is the complex space. We define as wj 2 R the j-th canonical basis vector of the complex space (a vector whose entries are all equal to 0 but 1 in the j-th position), j ¼ 1;   ; , representing the j-th complex. Let jR and jP be the indices of the reacting and the product complex of a given reaction, respectively (associated with the canonical basis vectors wjR and wjP , respectively). Such reaction can be described as the difference between wjP and wjR . We write:

k sout

Example: Consider the following biochemical reactions involving four compounds: C1 , C2 , C3 and C4 .

ð3bÞ ð3cÞ

ð6cÞ

ð8Þ

For the sake of clarity, we conclude the section by an elementary example.

3 C1

This decomposition is similar to the one used in [18]. The rank of N is denoted by: s ¼ rankðNÞ

N ¼ BS

ð3aÞ

In this way, we can decompose matrix N as N ¼ Nout  Nin :

Interestingly enough, the following relation holds:

Here three reactions take place. Namely, the reaction of association of three different molecules of C1 into C2 , the reaction of dissociation of species C2 into three molecules of C1 , and the reaction of association of species C2 and C3 into species C4 . Three moles of species C1 react to produce one mole of C2 at a rate of 1 1=½moles2 ½sec. Conversely, one mole of C2 dissociates into three moles of C1 at a rate of 2 1=½sec. Finally, one mole of C2 and two moles of C3 react to produce one mole of C4 at a rate of 3 1=½moles2 ½sec. The complexes involved in the network are collected in the vector X ¼ ½ X1 X2 X3 X4 0 , being X1 ¼ ð3 C1 Þ, X2 ¼ ðC2 Þ, X3 ¼ ðC2 ; 2 C3 Þ, X4 ¼ ðC4 Þ. We denote as c1 , c2 , c3 , c4 the concentrations of species C1 , C2 , C3 and C4 , respectively. We collect such variables in the vector c ¼ ½ c1 c2 c3 c4 0 . The four ( ¼ 4) complex vectors bj s are 4-dimensional (n ¼ 4), that is: b1 ¼ ½ 3 0 0 0 0 , b2 ¼ ½ 0 1 0 0 0 , b3 ¼ ½ 0 1 2 00 , b4 ¼ ½ 0 0 0 1 0 . The complex matrix is therefore: 2 3 3 0 0 0 60 1 1 07 7 B¼6 40 0 2 05 0 0 0 1 Now, according to their definition, we assign the vectors wj ’s (in the complex space), associated to the four complexes Xj ’s as w1 ¼ ½ 1 0 0 0 0 , w 3 ¼ ½ 0 0 1 0 0 , w4 ¼ w2 ¼ ½ 0 1 0 0 0 ,

581

An Observer for Mass-action Reaction Networks

½ 0 0 0 10 . We recall that the reacting complex (resp. product complex) of reaction 1 is X1 (resp. X2 ), the reacting complex (resp. product complex) of reaction 2 is X2 (resp. X1 ) and the reacting complex (resp. product complex) of reaction 3 is X3 (resp. X4 ). According to the definition (6) of vectors sk s (k ¼ 1; 2; 3), one obtains s1 ¼ w2  w1 , s2 ¼ w1  w2 , s3 ¼ w4  w3 . Matrix S in (7) is therefore: 2 3 1 1 0 6 1 1 0 7 7 S¼6 4 0 0 1 5 0 0 1

the exponential and the logarithmic operators are applied element-wise and ðÞ0 indicates transposition. Notice that, for mass-action networks, the number of elements of v is equal to the number of unknown model parameters, i.e., r ¼ p. In the sequel the following assumption is made.

According to (8) one can then compute the product B S and check that it coincides with the stoichiometric matrix 2 3 3 3 0 6 1 1 1 7 7 N¼6 4 0 0 2 5 0 0 1

where #B ðcÞ is defined as:

2.2. Rate Equations for Mass-action Networks Let cðtÞ 2 Rn indicate the vector whose entries are concentrations of chemicals (cðtÞ  0 element-wise). We define as xðtÞ 2 R the vector whose entries are the concentrations of the complexes of the network. The vector cðtÞ can be computed, in terms of xðtÞ, according to the following equation (see also [22]): cðtÞ ¼ B xðtÞ

ð9Þ

The rates of the reactions are indicated as vk ðk ¼ 1;    ; r). These rates depend upon the concentrations of the species cðtÞ 2 Rn and of the reaction rate constants  2 Rp . In general, the entries vk of such map are nonlinear functions vk ðc; Þ. The specific type of nonlinearity depends on the particular kinetic of each single reaction (e.g. one can consider the massaction law, Michaelis–Menten equation and Hill formula, see e.g. [8, 31]). Herein we focus on the mass-action kinetics, see e.g., [4, 22]. In such framework, the reaction rate is proportional to the concentrations of each reactant species. This leads to the relation [9]: 0

vðc; Þ ¼ diagðÞ cNin 0

where the symbol cNin has the following meaning   0 ~ cNin ¼ exp N0in logðcÞ

ð10Þ

Assumption 2: Each element of vector  is strictly positive. Notice that an alternative expression of formula (10) is, according to [22]: vðc; Þ ¼ diagðÞ S0in #B ðcÞ

0

#B ðcÞ ¼ cB 2 R

ð11Þ

ð12Þ

_ and xðtÞ _ Let cðtÞ represent the rate of change of cðtÞ and xðtÞ, respectively, at a given time t. Given the rate of the reactions vðÞ, they are computed as follows: x_ ¼ S vðc; Þ

ð13aÞ

c_ ¼ N vðc; Þ

ð13bÞ

We say that c is an equilibrium concentration if _ cðtÞ¼c ¼ 0. cðtÞj Correspondingly N vð c; Þ ¼ B S vð c; Þ ¼ 0. Of course, this does not necessarily imply that S vð c; Þ ¼ 0 (leading to x_ ¼ 0). When, in a given equilibrium point c, also x_ ¼ 0 holds true, we say that c is complex balanced. Biochemical reaction networks are in general stimulated by exogenous signals. They can be imposed through experiments, generated by other cells or inside the cell itself. We consider such quantities as inputs for the system, provided that their values can be measured. In this paper we consider inputs whose effect enter the model ‘‘linearly’’. Consistent with biochemical practice, we assume input values to be non-negative. In in-vitro experiments, a single concentration or a single reaction rate is often measured. In the sequel we will consider as outputs of the model single species concentrations or single reaction rates. In conclusion, in view of (13b) and (10), our model can be written as: 0

c_ ¼ N diagðÞ cNin þ G u  y¼

yc yv



 ¼

Hc 0

0 Hv

  c v

ð14aÞ ð14bÞ

The vector u is a pu -dimension vector, and matrix G 2 Rnpu is a matrix describing the effect of the inputs on

582

M. Farina and S. Bittanti

the rate of variation of single species concentrations. In the simplest case the n-dimension columns of matrix G are canonical basis vectors of Rn . Further exogenous influence terms could be introduced, other than linear terms [3]. According to the assumption made, the rows of the matrices Hv and Hc are vectors with only one non-zero element. We denote with pc and pv the number of measured species and reaction rates, respectively. Note that, in some experimental setups, the sum of several species is measured, since compounds can be involved in many complexes. The extension of the results of this paper to this case will be subject of further investigation.

introduce a further assumption concerning the possible boundary equilibria (see [4]). We define the stoichiometric space as D ¼ _ ¼ spanfn1 ;    ; nr g 2 Rn . Intuitively, it can be spanfcg thought as the set representing the constraints to which all the states must attain, since all variations in c must happen according to the proportions dictated by the reaction vectors. The scalar s, defined as the rank of the stoichiometric matrix N, (5), is the dimension of such subspace. Finally we define, corresponding to each initial condition c0 2 Rn0 , the sets: S c0 :fc0 þ d=d 2 Dg \ Rn0 is called reaction simplex n Sþ c0 :fc0 þ d=d 2 Dg \ R>0

2.3. Properties of Equilibria of Mass-action Chemical Networks

is called positive reaction simplex

In this section we recall one of the main results on the analysis of mass-action networks, the celebrated zero deficiency theorem. This theorem deals with the equilibria of a sub-class of such networks, specifically the socalled zero deficiency networks. In order to introduce this class, some preliminary concepts are presented, such as the notions of linkage class, weak reversibility and deficiency. The reader who is already familiar with such concepts and results can skip this section. Directed graphs are suitable graphical tools for reaction networks. The set of indices J ¼ f1;    ; g (related to complexes) can be subdivided into subsets Li (linkage classes) of connected components. We say that the complex j1 and complex j2 belong to the same linkage class (j1  j2 ) if there exists a path (not necessarily a directed path) which connects the complexes j1 to j2 . Definition 1: (Weak Reversibility) A network is said to be weakly reversible if there exist a directed path between each element of a linkage class and any other element of the same linkage class.

ð16aÞ

ð16bÞ

Here we denote Rn0 the positive orthant, and Rn>0 the strictly positive orthant of the state-space. One can show that S c0 and S þ c0 are forward invariant. i.e., if cð0Þ 2 S c0 ; cðtÞ 2 S c0 8 t  0. The subspace S c0 can be thought as the set of all the attainable points, starting at c0 as initial condition. It is also referred to, in the literature, as stoichiometric compatibility class. Assumption 3: The system (13b) has no equilibria on the boundary of the positive reaction simplex S þ c0 . Theorem 1 (Zero Deficiency Theorem): For a zero deficiency and weak reversible network satisfying Assumption 1, – there is a unique equilibrium c in the positive reaction simplex S þ c0 ; – the equilibrium c is asymptotically stable; – the equilibrium c is globally asymptotically stable if and only if Assumption 3 holds. Otherwise, the equilibrium c is only locally asymptotically stable. &

Definition 2: The deficiency of a network is defined as:  ¼ls

ð15Þ

where  is the number of complexes, l is the number of linkage classes and s is defined as in (5). It can be shown that   0, see [22]. The celebrated zero deficiency theorem deals with the existence, uniqueness and stability of the equilibria in zero deficiency networks. Such result has been stated in a number of different forms, starting from the 1970s, see [12, 14, 23–25]. Here we recall the most comprehensive statement [34]. To this purpose we

2.4. Degradation and Formation Reactions Formation and degradation reactions are widely studied in the literature. Nevertheless they are usually neglected in the studies on mass-action network models because they may induce the violation of some main assumptions (e.g. weak reversibility). Here we include them in our modelization. The peculiarity of such reactions is that either the reactant (XR ) or the product (XP ) complex is not explicitly present. That is, there exist three main types of reactions, i.e.:

583

An Observer for Mass-action Reaction Networks

XR ! XP

\classical" reaction

   ! XP XR !   

formation reaction degradation reaction

Degradation. Consider a complex (denoted with index j) consisting of a single species ( i). Reaction k is  said to be a degradation reaction if nkout ¼ 0. That is, the product species lacks. From this follows that i-th canonical nk ¼ bj . Clearly, bj ¼ ei being ei the  basis vector in Rn . The dynamics of a network may vary significantly (from a qualitative point of view) by adding degradation terms. Consider for example a zero deficiency weakly reversible network. If degradation terms are neglected, a steady state equilibrium in the positive state space is guaranteed by Theorem 1. When degradation is included in the model, the deficiency may continue to be equal to zero1 but the property of weak reversibility is necessarily lost2and a strictly positive steady state equilibrium in the positive state space is not guaranteed to exist (see Theorem 1). Formation. Consider a complex (denoted with index j) consisting of a single species ( i) forming with a rate uj .  Reaction k is said to be a formation reaction if nkin ¼ 0. That is, the reactant species lacks. From this follows i-th canonical that nk ¼ bj . Clearly, bj ¼ ei being ei the  basis vector in Rn . In this framework one can consider the formation rates uj as reaction rate (additional element of vector v, independent of the variables c). If the value of uj is known, and since it assumes positive values, an equivalent formulation is to consider such terms as inputs. That is, being uj the rate of formation of the species ci , we can write, recalling (13b): c_i ¼ e0i N vðc; Þ þ uj where ei is the i-th canonical basis vector (with absolute value equal to 1) of Rn . Remark that, in this case,

1

2

For the computation of the deficiency , one should consider a ‘‘virtual product complex’’ j0 associated to such a reaction. So the total number of complexes tot increases, as well as the number of reactions involved and the dimensions (e.g. the rank s) of matrix N. Note that also the number of linkage classes may increase. This implies that, in a network, if degradation terms are considered, the value of  may be invariant with respect to a case when degradation is neglected. Indeed, if we represent degradation as an irreversible reaction occurring between a ‘‘real’’ complex j and a ‘‘virtual’’ one j0 , there is no directed path from j0 to any other complex of the linkage class. Hence, weak reversibility is not satisfied.

formation is not regarded as a reaction, and the total number of reactions r is conserved.

3. Main Result: The Observer We now tackle the problem of estimating the concentrations of species in a mass-action network from the available measurements. For model (14), it is apparent that this is just an observability problem, namely state estimation from output measurements. First we consider the problem for a weakly reversible mass-action network (Section 3.1). In Section 3.2, we specialize to zero deficiency networks. Finally, in Section 3.3, we encompass formation and degradation effects. 3.1. The Observer for a Weakly Reversible Mass Action Network The problem tackled here is twofold. When the network and the experimental measurement facilities are given, the basic issue is how to process the available data in order to estimate the concentrations with high accuracy, if possible. This means designing a suitable data processor (observer synthesis). The problem can also be posed the other way round, namely to select the measured variables so as to maximize the information contained in data for the most accurate estimation of the concentrations (experimental design). In this paper, both the problem of experimental design and that of observer synthesis are tackled, given the model (14). As for the latter, the main result is that, under a given condition, we guarantee that the concentrations estimated by means of our Luenberger observer are asymptotically correct, namely they tend to the true values. The condition is a limitation on the rates of variation of the concentrations of the complexes (i.e. the variable xðtÞ). Note that, such condition is automatically satisfied for systems of zero deficiency weak reversible class, in steady state. In order to derive the observer equations, we first re-write system (14): 8 0 < c_ ¼ N diagðÞcNin þ G u y ¼ Hc c : c 0 yv ¼ Hv diagðÞcNin Taking advantage of the system structure, we consider the filter: 0 c^_ ¼ N diagðÞ c^Nin þ G u þ Kc ðyc  Hc c^Þ  0 þ Kv yv  Hv diagðÞ c^Nin

ð17Þ

584

M. Farina and S. Bittanti

as observer. Here c^ðtÞ is the estimate of concentration cðtÞ. The design of the observer amounts to selecting the gain matrices Kc and Kv so as to guarantee the asymptotic convergence of c^ðtÞ to cðtÞ at least locally. To be precise, we define the state estimation error eðtÞ as: eðtÞ ¼ c^ðtÞ  cðtÞ

ð18Þ

We say that the observer is locally asymptotically convergent if, for each initial condition for the observer lying in a neighborhood of the initial true concentration, eðtÞ ! 0. Theorem 2: For a weakly reversible mass-action network, if the following two conditions hold:

  B ¼n ð19aÞ rank Hc

 rank

S0 Hv S0in

 ¼

ð19bÞ

then, by taking: Kc ¼ H0c

ð20aÞ

Kv ¼ Nout H0v

ð20bÞ

the observer (17) converges locally for sufficiently low rates of variation of the complex concentrations xðtÞ _ (namely under condition kxðtÞk < x for a suitable  value of x ) for strictly positive values of species concentrations (cðtÞ > 0 element-wise). & Since matrices B, S and Sin are determined by the network topology, the possibility of satisfying the observability conditions (19) entirely depends on the matrices Hc and Hv (that is on the selected species concentrations and reaction rates to be measured). In this regard, it is worthy considering that a main issue is the minimization of the number of outputs. Both the constraints on the initial values of the observer state (that is the bound imposed on e(0)) and on the derivative of the complex concentrations xðtÞ are made explicit in the proof in Appendix A. Note that, _ under certain assumptions, the condition kxðtÞk < x can be replaced by an assumption on the dynamics of _ the species concentrations (kcðtÞk < c for a suitable c ). This is discussed in the following paragraphs. Finally, note that, since Theorem 2 holds if the considered network is weakly reversible, it does not encompass networks where species are affected by degradation which is, by definition, a non reversible reaction. This case is studied in Section 3.3.

3.2. The Case of Zero-deficiency Networks For weakly reversible zero deficiency networks, the condition on the dynamics of variable xðtÞ can be expressed as conditions on the dynamics of cðtÞ, as in the following corollaries of Theorem 2. Corollary 1: Consider a zero deficiency weakly reversible network. Under the assumptions (19) and with choice (20) of the observer gains, the observer (17) converges locally for sufficiently low rates of variation of strictly positive values of the species concentrations _ cðtÞ (namely under condition kcðtÞk < c for a suitable value of c and cðtÞ > 0 element-wise). Finally, a further corollary can be stated. The case _ ¼ 0) for treated here is the steady state (that is cðtÞ zero deficiency weakly reversible networks (with u ¼ 0). Corollary 2: Consider a zero deficiency weakly reversible network in an equilibrium condition with no exogenous signals (i.e. u ¼ 0). Assume Assumption 3 be satisfied. Under the assumptions (19) and with choice (20) of the observer gains, the observer (17) converges locally.

3.3. Degradation and formation reactions As discussed in Section 2.4, formation and degradation are particular reactions where either the reactant complex or the product complex is not specified. In this section we show that the observer designed disregarding degradation or formation is also suitable for state estimation when such phenomena occur, provided that the rates of formation are known. Corollary 3: Given a system (14), where the reaction scheme obtained disregarding degradation and formation of species is weakly reversible, the topology of which is described by matrices B0 , S0in and S0 . If Hc and Hv are computed according to (19) (with B0 , S0in and S0 ), then the observer (17) converges locally even if the species are affected by formation and degradation reactions, for sufficiently low rates of variation of the complex concentrations xðtÞ (namely under condi_ tion kxðtÞk < x for a suitable value of x ) for strictly positive values of species concentrations (cðtÞ > 0 element-wise). The proof of the statement above follows from the proof of Theorem 2, and is given in Appendix C.

585

An Observer for Mass-action Reaction Networks

3.4. Discussion Our main result provides an observer which can be applied to a fairly general class of dynamical biochemical network models, the weakly reversible massaction systems. First we have established a formalism for a handy description of such class of models, which relies on previous works on biological systems and enzyme kinetics modelling [5, 13–15, 17, 22]. Thanks to such a formalism, we can provide some conditions of application of a Luenberger observer of the type (17). To be precise, with a Lyapunov-type argument, we can show that the selection of gains (20) not only guarantees convergence of the state estimate but also supplies observability conditions which are independent from the equilibrium point. Thus we can prove a tool for experiment design of general and easy applicability. Once the observability conditions are established in this way, both different choices of the observer gains and alternative observers could be employed, in principle, to address specifications, such as increase of estimation speed and better noise rejection. Future work will be devoted to the assessment of the performance achievable by applying different techniques for the observer synthesis, ranging from extended Kalman filters to other nonlinear system observers, [19]. An extension of the latter is described in [6, 20], where the bottleneck of transforming the system in canonical form is avoided and more flexibility for the choice of the observer gains is provided. In Section 3.2, the results of Section 3.1 are applied to a subclass of mass-action networks, the so-called zero-deficiency networks. This class is explored also in the paper [4], where conditions for detectability are given and an observer is worked out, provided that such conditions are satisfied. In [4], an equivalent formulation of detectability is provided accounting for the steady state of the system.

4. A Case Study: the Edelstein Network In this section we apply the observer worked out in this paper to a ‘‘benchmark’’ reaction network, previously referred to in [11] as the Edelstein network, originally investigated in [7]. This reaction network is of great interest, since the induced differential equations admit multiple positive equilibria for some values of the rate constants. In particular, we now take into account of the case discussed in [11]. To be precise, here we study: (a) the Edelstein network (Section 4.2); (b) the ‘‘modified’’ Edelstein network presented in Problem 3.A.2 and Example 3.2.3 in [11], respect-

ively (Section 4.2); (c) a version of the plain Edelstein network, affected by species degradation (Section 4.3)3. For simplicity of notation, the names of the reacting species have been changed with respect to [11]

4.1. Plain Edelstein network The reaction scheme is the following: k1

C1 Ð 2C1 k1

k3

k2

C1 þ C2 Ð

C3 Ð C2

k2

k3

The state variables are the n ¼ 3 concentrations of C1 , C2 and C3 , c ¼ ½ c1 c2 c3 0 . The reaction network presents r ¼ 6 reactions (with constant reaction rates k1 , k1 , k2 , k2 , k3 and k3 ). Such reaction rates are collected in the vector v ¼ ½ k1 c1 ; k1 c21 ; k2 c1 c2 ; k2 c3 ; ¼ ½ v1 v2 v3 v4 v5 v6 0 :

0

k3 c3 ; k3 c2 

As suggested by [11], we consider the set of parameters  ¼ ½8:5; 1; 1; 1; 1; 0:2, such that the system admits multiple positive equilibria. Note that the reaction network is weakly reversible. Matrices B and S are: 2

1 B ¼ 40 0

2 1 0 1 0 0

0 0 1

2

1 6 1 0 6 15; S ¼ 6 6 0 4 0 0 0 3

1 1 0 0 0

0 0 0 0 1 1 1 1 0 0

0 0 0 1 1

3 0 0 7 7 0 7 7 1 5 1

The network presents  ¼ 5 complexes involved in the reactions, specifically X1 ¼ ðC1 Þ, X2 ¼ ð2C1 Þ, X3 ¼ ðC1 ; C2 Þ, X4 ¼ ðC3 Þ, X5 ¼ ðC2 Þ. Since  ¼ 5, l ¼ 2 and s ¼ rankðB SÞ ¼ 2, the deficiency of the network is  ¼ 5  2  2 ¼ 1. To satisfy conditions (19a) and (19b), the output measurement matrix Hc can be assigned as an empty matrix, while Hv can be assigned, for example, as   1 0 0 0 0 0 : Hv ¼ 0 0 1 0 0 0 3

In the following examples the initial conditions cð0Þ (for the state variables) and c^ð0Þ (for the estimates) are assigned arbitrarily: cð0Þ ¼ ½ 12:3; 9:6; 16:5  c^ð0Þ ¼ ½ 0:001; 0:001; 0:001 

586

M. Farina and S. Bittanti

4.2. Modified Edelstein network In Example 3.2.3 of [11], a modified version of the Edelstein network is presented, whose state variables are proved to have only one equilibrium point (see the interesting discussion on the deficiency one theorem in [11]). The reacting species and complexes are the same as in Section 4.1, but the number of reactions is now larger: r ¼ 8 in place of r ¼ 6 as in the plain network. So, matrix B remains unchanged, while matrix S is now: 2 3 1 1 0 0 0 0 1 1 6 1 1 0 0 0 0 0 0 7 6 7 6 S¼6 0 0 1 1 0 0 0 0 7 7 4 0 0 1 1 1 1 1 1 5 0 0 0 0 1 1 0 0 We assign k4 ¼ 0:2 and k4 ¼ 1. As in the previous case, to satisfy (19a) no species concentration measurements are required, while the number of measured reactions should now be one, to satisfy (19b). In fact, Hv ¼ ½1 0 0 0 0 0 0 0 satisfies (19b), and the vector of measurements consists in yv ¼ v1 solely. Fig. 2 shows the trajectories of the species concentrations and their estimates.

4.3. Plain Edelstein Network Affected by Species Degradation As an application of Corollary 3, we now assume that degradation affects species C3 , with degradation rate kd3 ¼ 1. As stated in Corollary 3, in the case where degradation does not affect the network species, the measurements consist in y ¼ ½ v1 v3 0 . Fig. 3 shows the convergence of the proposed observer. Note that, despite the system’s steady state is on the boundary of the stoichiometric compatibility class, convergence is guaranteed as long as c > 0 and the bound of the convergence rate are verified. Nonetheless, the verification of the latter has not been explicitly addressed.

5. Conclusions Huge is the variety and complexity of phenomena studied in life sciences. In the last decades, mathematical modelling methods have provided useful tools for the comprehension of some of them, in particular, within systems biology, the dynamics of 20 18 Species concentrations (arbitrary units)

This is equivalent to say that no species measurements are required, and the measurement output is simply given by: y ¼ ½ v1 v3 0 . When the observer developed in Section 3.1 is used, the obtained estimates states are depicted in Fig. 1, where the plot of the true state diagram is also given for ease of comparison with the estimates.

16 14 12 10 c (t) 1

8 6 4

c (t) 2

c3(t)

2 0 0

2

4

6

Time (arbitrary units)

8

10

Fig. 2. State trajectory of the system (solid lines) and of the observer (dashed lines) for the modified Edelstein network. 20

18

16 Species concentrations (arbitrary units)

Species concentrations (arbitrary units)

20

c (t) 3

18

14 12 10 c2(t)

8 6

c (t)

4

1

2 0 0

1

2

3

4

5

Time (arbitrary units)

Fig. 1. Evolution of the system state variables (solid lines) VS their estimates (dashed lines) obtained by the observer developed in Section 3.1 for the plain Edelstein network.

16 14 12 10

c1(t)

8 6 4

c (t) 3 c (t) 2

2 0 0

1

2

3

4

5

Time (arbitrary units)

Fig. 3. State trajectory of the system (solid lines) and of the observer (dashed lines) in the example of Section 4.3.

An Observer for Mass-action Reaction Networks

cellular networks. In this framework a central problem is the determination of species concentrations in biochemical reaction networks. The question of estimating the (unmeasurable) state of a system from measured signals is known in control science as the observer problem. In this paper we exploit control science methods to work out a nonlinear observer providing fair estimates of the concentrations for a class of mass-action biochemical networks. Under certain assumptions, we prove that the state estimates provided by the observer converge to the true concentrations as the number of available data tends to infinity. Finally we discuss the application to zero deficiency mass-action networks, in particular to the so-called Edelstein network as a final case study. The observer methodology can also be used to address the question of experiment design, namely to envisage the set of variables which must be measured to guarantee the possibility of estimanting the unknowns (observability issue).

Acknowledgments Research has been supported by The Italian National Research Project ‘‘New Techniques of Identification and Adaptive Control for Industrial Systems’’, and by CNR-IEIIT. The authors would like to thank the Reviewers for their constructive comments and criticism.

References 1. Bastin G, Dochain D. Online Estimation and Adaptive Control of Bioreactors. Elsevier, 1990 2. Bastin G, van Impe JF. Nonlinear and adaptive control in biotechnology: a tutorial. Eur J Control 1995; 1: 1–37 3. Chaves M, Sontag ED. State-estimators for chemical reaction networks of Feinberg-Horn-Jackson type. Eur J Control 2002; 8: 343–359 4. Chaves M. Input to state stability of reaction-controlled chemical networks. SIAM J Control Optim 2005; 44(2): 704–727 5. Cornish-Bowden A. Fundamentals of Enzyme Kinetics, 3rd edn. Portland Press, 2004 6. Dalla Mora M, Germani A, Manes C. Design of state observers from a drift-observability property. IEEE Trans Autom Control 2000; 45(8): 1536–1540 7. Edelstein B. Biochemical model with multiple steady states and hysteresis. J Theor Biol 1970; 29(1): 5762 8. Fall CP, Marland ES, Wagner JM, Tyson JJ. Computational Cell Biology. Springer, 2002 9. Farina M, Findeisen R, Bullinger E, Bittanti S, Allgo¨wer F, Wellstead P. Results towards identifiability properties of biochemical reaction networks. Proceedings of Conference on Decision and Control, San Diego, CA, USA, 2006, pp. 2104–2109

587 10. Farina M, Findeisen R, Bullinger E, Bittanti S. An observability based strategy for parameter identification in systems biology. Foundations of System in Engineering (FOSBE) Conference 2007, Stuttgart, DE, 9–12/09/ 2007, 2007 11. Feinberg M, Horn FJM. Dynamics of open chemical systems and the algebraic structure of the underlying reaction networks. Chem Eng Sci 1974; 29: 775 –787 12. Feinberg M. Complex balancing in general kinetic systems. Arch Ration Mech Anal 1972; 49: 187–194 13. Feinberg M. Mathematical aspects of mass action kinetics. In: Lapidus L, Amundson N (eds), Chemical Reaction Theory: A Review. Prentice Hall, NJ, 1977 14. Feinberg M. Chemical reaction network structure and the stability of complex isothermal reactors I. The deficiency zero and the deficiency one theorems. Chem Eng Sci 1987; 42: 2229–2268 15. Feinberg M. Chemical reaction network structure and the stability of complex isothermal reactors ii: multiple steady states for networks of deficiency one. Chem Eng Sci 1988; 43: 1–25 16. Feinberg M. The existence and uniqueness of steadystates for a class of chemical reaction networks. Arch Ration Mech Anal 1995; 132: 311–370 17. Feinberg M. Lectures on chemical reaction systems. Lectures delivered at the Mathematics Research Center, University of Winsconsin-Madison in the autumn of 1979. Website: http://www.che.eng.ohio-state.edu/ FEINBERG/LecturesOnReactionNetworks/ 18. Gatermann K, Huber B. A family of sparse polynomial systems arising in chemical reaction systems. J Symbol Comput 2002; 33: 275–305 19. Gauthier JP, Hammouri H, Othman S. A simple observer for nonlinear systems. applications to bioreactors. IEEE Trans Autom Control 1992; 37(6): 875–880 20. Germani A, Manes C. State observers for nonlinear systems with smooth/bounded input. Kybernetika 2000; 36(1): 31–42 21. Gilles ED. Control—key to better understanding biological systems. Automatisierungstechnik 2002; 50: 7–17 22. Gubermann JM. Mass action reaction networks and the zero deficiency theorem. Master’s thesis, Harvard University, Cambridge, 2003 23. Horn FJM, Jackson R. Stability and complex balancing in mass-action systems with three short complexes. Arch Ration Mech Anal 1972; B47: 81–116 24. Horn FJM. A connection between stability and graphs in chemical kinetics – II: stability and the complex graph. Proc R Soc 1973; A334: 313–330 25. Horn FJM. A connection between stability and graphs in chemical kinetics – I: stability and the reaction diagram. Proc R Soc 1973; A334: 299–312 26. Horn FJM. Stability and complex balancing in massaction systems with three short complexes. Proc R Soc 1973; A334: 331–342 27. Khalil H. Nonlinear systems, 3rd edn. Prentice Hall 28. Kholodenko B, Sauro HM, Westerhoff HV. Control by enzymes, coenzymes and conserved moieties. A generalization of connectivity theorem of metabolic control analysis. Eur J Biochem 1994; 125: 179–186 29. Kitano H. Looking beyond the details: a rise in systemoriented approaches in genetic and molecular biology. Curr Genet 2002; 41: 1–10

588

M. Farina and S. Bittanti

30. Kitano H. Systems biology: a brief overview. Science 2002; 295: 1662–1664 31. Klipp E, Herwig R, Kowald A, Wierling C, Lehrach H. Systems Biology in Practice: Concepts, Implementation and Application. Wiley-VCH, 2005 32. Pous NM, Rajab A, Flaus A, Engasser JM, Cheruy A. Comparison of estimation methods in biotechnological processes. Chem Eng Sci 1914; 43: 1909–1914 33. Sauro HM. Structural Analysis, Part II: Conserved Relations. Study notes, Version 1.0, January 2000, 2000 34. Sontag ED. Molecular systems biology and control. Eur J Control 2005; 11: 396–435 35. Sontag E. Structure and stability of certain chemical networks and applications to the kinetic proofreading model of T-cell reaceptor signal transduction. IEEE Trans Autom Control 2001; 46(7): 1028–1047 36. Wellstead P. The dynamic systems approach to control and regulation of intracellular networks. FEBS Lett 2005; 579: 1846–1853 37. Wellstead P. A plea for more theory in molecular biology. In: Butcher E, Bringmann P, Weiss B (eds), Systems Biology – Applications and Perspective, SpringerVerlag, 2006 38. Wellstead P. Control opportunities in systems biology. Plenary Address, IFAC Symposia, CAB 2007 and DYCOPS 2007 39. Wolkenhauer O, Cho KH, Kolch W. Systems biology: towards and understanding of the dynamics of life. Bioforum Eur 2004; 1: 50–52

Appendix A A. Proof of the main Theorem 2 In this appendix we prove the main statement of the paper, namely Theorem 2, given in Section 3. The proof requires a number of preliminary results, and it is therefore organized in a number of steps. Basically we resort to the Lyapunov theory and we provide conditions which guarantee the asymptotic convergence of the observer. To achieve this objective we require that a Lyapunov-type matrix inequality be satisfied and that the initial conditions of the observer state lie in a neighborhood of the initial true concentration (local result). Such inequality holds true provided that the rank conditions (19) on the output measurement matrices Hc and Hv be satisfied. A further constraint on the rates of growth of the complex concentration xðtÞ arises from the study. First of all some notation is introduced. A.1. Notation We define the following matrices:  0  N A1 ¼  in 2 Rðrþpc Þn Hc

 A2 ¼

ðN  Kv Hv Þ0 K0c



2 Rðrþpc Þn

ð21bÞ

D1 ðcÞ ¼ diagðcA1 Þ 2 Rðrþpc Þðrþpc Þ

ð21cÞ

D2 ðcÞ ¼ ðdiagðcÞÞ1 2 Rnn

ð21dÞ

 D3 ¼

diagðÞ 0pc r

0rpc Ip c



2 Rðrþpc Þðrþpc Þ

ð21eÞ

AðcÞ ¼ A02 D3 D1 ðcÞ A1 D2 ðcÞ 2 Rnn

ð21fÞ

where Ik indicates the identity matrix of dimension k, and 0ij indicates a ði  jÞ matrix whose entries are all equal to zero. A.2. A Lyapunov Matrix Inequality Guarantees Local Convergence. Recalling that cðtÞ denotes the state of the system and that c^ðtÞ denotes its estimates provided by the observer, the estimation error is eðtÞ ¼ c^ðtÞ  cðtÞ, as already stated in equation (18). The observer is said to be locally convergent if there exists a set D such that 0 2 D and keðtÞk ! 0 when eðt0 Þ 2 D at a certain time instant t0 . In this section we prove that the local convergence is ensured when a Lyapunov-type inequality holds true. To this purpose, we first study the dynamics of eðtÞ. Dynamics of the estimation error eðtÞ. From (14) and (17): e_ ¼ ðN  Kv Hv Þ ðvðcðtÞ þ e; Þ  vðcðtÞ; ÞÞ  Kc Hc e ¼ ½ðN  Kv Hv Þ diagðÞ 0  0   0 1 Nin Nin B Hc  cðtÞ Hc C C  Kc B @ðcðtÞ þ eÞ A ð22Þ The last term in (22) will be denoted as "ðÞ. Recalling definitions (21) we write:  "ðe; cðtÞÞ ¼ ðcðtÞ þ eÞA1  cðtÞA1 By the Taylor series:

ð21aÞ

"ðe; cðtÞÞ ¼  diagðcA1 Þ A1 ðdiagðcÞÞ1 e þ rðe; cÞ

ð23Þ

589

An Observer for Mass-action Reaction Networks

where rðe; cÞ is the Peano remainder. Remark that rðe; cÞ is a continuous function of e such that: lim

e!0

krðe; cÞk ¼0 kek

ð24Þ

Eq. (22) can be rewritten as: e_ ¼ A02 D3 "ðe; cðtÞÞ

ð25aÞ

¼ AðcÞ e þ A02 D3 rðe; cÞ

ð25bÞ

In eq. (25b)) one can distinguish the linear parameterically varying term AðcÞ e. Here c plays the role of the varying parameter. Intuitively, around the origin, the behaviour of the system is lead by the linear part. We are now in a position to prove the following result. Lemma 1: The observer (17) is locally convergent if there exists a symmetric positive definite matrix (whose value may depend on cðtÞ) QðcÞ 2 Rnn such that the following inequality holds4: _ QðcÞ þ AðcÞ QðcÞ þ QðcÞ AðcÞ0 > 0

ð26Þ &

Proof: Let the following quadratic function be the candidate Lyapunov function: 1 Vðe; cðtÞÞ ¼ e0 QðcÞ1 e 2 For simplicity, we denote PðcÞ as the inverse of QðcÞ, i.e., PðcÞ ¼ QðcÞ1 . Since QðcÞ is a symmetric definite positive parameter-dependent matrix by definition, VðÞ is positive definite. The derivative of Vðe; cðtÞÞ is computed as: _ e þ e0 PðcÞe_ _ cðtÞÞ ¼ 1 e0 PðcÞ Vðe; 2 1 _  PðcÞ AðcÞ  AðcÞ0 PðcÞ e ¼ e0 PðcÞ 2 þ e0 PðcÞA02 D3 rðe; cÞ Assume that there exists PðcÞ such that the matrix _ SðcÞ ¼ PðcÞ  PðcÞ AðcÞ  AðcÞ0 PðcÞ is negative definite. If this holds, the term 12 e0 SðcÞ e is negative for any value of eðtÞ not identically equal to zero. Then V_ is negative if at least one of the following conditions is satisfied: 4

_ The matrix QðcÞ is defined as follows. Being qij ðcÞ the entries of @ _ qij ðcÞ @c QðcÞ, the entries of QðcÞ are @c @t

IÞ e 0 PðcÞA02 D3 rðe; cÞ 0 IIÞ e0 PðcÞA02 D3 rðe; cÞ < 12 e0 SðcÞ e Now, by definition: 0 – e PðcÞA02 D3 rðe; cÞ kek2 krðe; cÞk2 kPðcÞA02 D 3 k2 – 12 e0 SðcÞ e  12 kek22 jmin ðSðcÞÞj, being min ðSðcÞÞ the eigenvalue with minimum magnitude of SðcÞ. Property II) is then guaranteed if: krðe; cÞk2 1 jmin ðSðcÞÞj < kek2 2 kPðcÞA02 D3 k2

ð27Þ

Since rðe; cÞ is a continuous function of e and (24) holds, there exists a value e > 0 (function of variable c) such that condition (27) is satisfied for all e such that kek2 < e. The value e defines a region around the origin such that V_ < 0, provided that the negative definiteness of SðcÞ is ensured. This guarantees the existence of a region of attraction D such that 0 2 D and eðtÞ ! 0 if, at initial time instant t0 , eðt0 Þ 2 D. Recall that QðcÞ ¼ PðcÞ1 and note that _ _ QðcÞPðcÞQðcÞ ¼ QðcÞ. We come then to the conclusion that there exists PðcÞ > 0 such that SðcÞ < 0 iff there exists QðcÞ > 0 such that condition (26) holds. &

A.3. The Observer Gains Kc and Kv In this section we analyze the inequality (26). Here we take the observer gains Kc and Kv as in (20). Under specific assumptions on Hc and Hv we can guarantee (26), specifically the rank conditions (19). In the sequel we will use the notion of quasi-positive definite matrix. Recall that L is said to be a quasi-positive definite if L þ L0 is positive definite. In the following lemma we make the assumption that Sv 0. Notice that such condition is satisfied by the state variables in a stationary condition, both by an autonomous system (i.e., u ¼ 0, where Sv ¼ 0) and by a system endowed with non negative input variables (i.e., since u  0, if x_ ¼ 0 then Sv 0). Lemma 2: For a weakly reversible network, let the rank conditions (19) be satisfied and consider two vectors v 2 Rr and yc 2 Rpc strictly positive element-wise. Assume that the Luenberger gains be as in (20). Let Sv 0

ð28Þ

590

M. Farina and S. Bittanti

element-wise. Then matrix:  0 diagðvÞ Lðv; yc Þ ¼ A2 0pc r

and Sout have, by construction, only one nonzero element. So:



0rpc A1 diagðyc Þ

is quasi-positive definite.

ð29Þ &

Proof: According to the definition of a quasi definite positive matrix Lðv; yc Þ, in (29), is quasi definite positive iff, for all vectors  2 Rn not identically equal to zero, the following holds: Jð; v; yc Þ ¼ 0 Lðv; yc Þ  > 0  v (the ‘‘orthogonal matrix’’ of Hv ) be Let H defined such that both the following properties are satisfied:  v are canonical basis vectors

the row vectors of H of Rr  v are the vectors of the basis of

the row vectors of H kerðHv Þ

diagðS0in &Þ S0in & ¼ S0in diagð&Þ &;

ð31aÞ

diagðS0out &Þ S0out & ¼ S0out diagð&Þ &

ð31bÞ

Note that, if the network is weakly reversible, v 6¼ 0 element-wise satisfying Sv ¼ 0 exists. Thank to Assumption (28), it holds that: v0 S0 0 ) v0 S0in diagð&Þ &  v0 S0out diagð&Þ & ð32Þ Hence, (13) and (32) imply that:   0 H  0 Jv ð&Þ ¼ v0 diagðS0in &Þ Sin  Sout H v v & 1  v0 fð&Þ 2

ð33Þ

where: Notice that, since we assume the rows of Hv to be canonical basis vectors of Rr , it is always possible to  v satisfying these conditions. We compute: find H Jð; v; yc Þ ¼ 0 Lðv; yc Þ  n  0 0 H  ¼ 0 Nin  Nout H v v diagðvÞ Nin o þ H0c diagðyc Þ Hc  n 0 H  ¼ 0 H0c diagðyc Þ Hc  þ 0 Nin  Nout H v v o diagðvÞ N0in  We denote: Jv ðÞ ¼ 0



  0 0 H  Nin  Nout H v v diagðvÞ Nin  ð30aÞ

Jc ðÞ ¼ 0 H0c diagðyc Þ Hc

ð30bÞ

h   0 H  0 fð&Þ ¼ diagðS0in &Þ Sin  Sout H v v &   i 0 H v 0& þ diagðS0 &Þ Sout  Sin H out

v

The vector fð&Þ has r elements. Given the k-th reaction, we denote rk and pk the indices of the reacting complex and of the product complex, respectively. According to this notation, and (33), the k-th element of fð&Þ is: &r2k þ &p2k  2 k &rk &pk ;

ð34Þ

where k ¼ 0 if ek (the k-th canonical basis vector) is a row of matrix Hv and k ¼ 1 otherwise. Notice that the elements (34) are non negative. Since v > 0 element-wise (since c > 0 element-wise), Jv ðÞ  0. Note that, by construction, Jc ðÞ  0. We first conclude that: Jð; v; yc Þ  0

ð35Þ

Notice that:   0 0 H  Jv ðÞ ¼ 0 Nin  Nout H v v diagðNin Þ v: We also denote & ¼ B0 . In this way, by recalling (8), we obtain:   0 H  0 Jv ðÞ ¼ v0 diagðS0in &Þ Sin  Sout H v v & Notice that the entries of matrix Sin (respectively Sout ) are either 1 or 0. Moreover, the columns of Sin

In the sequel, to conclude the proof, we show that Jð; v; yc Þ ¼ 0 implies that  ¼ 0. First, Jð; v; yc Þ ¼ 0 only if &r2k þ &p2k  2 k &rk &pk ¼ 0 for all k ¼ 1;    ; r

ð36Þ

Equation (36) is equivalent to &rk ¼ &pk for all reactions k ¼ 1; ::: ; r and, furthermore, that &rk ¼ 0 for all measured reactions. In other terms:

591

An Observer for Mass-action Reaction Networks



 S0 &¼0 Hv S0in

ð37Þ

If Assumption (19b) is verified, equation (37) holds only if & ¼ 0. So, &rk ¼ 0 for all k ¼ 1;    ; r, which corresponds to B0  ¼ 0

ð38aÞ

On the other hand, Jc ¼ 0 implies that Hc  ¼ 0

ð38bÞ

Now it is possible to prove the main Theorem 2. Proof (Theorem 2): From Lemma 1, the observer (17) is locally asymptotically convergent if there exists QðcÞ 2 Rnn > 0 such that the inequality (26) holds. Take QðcÞ ¼ ðD2 ðcÞÞ1 ¼ diagðcÞ, then the following inequality to be satisfied: ð39Þ

Take Hv such that (19b) is satisfied. For an autonomous system (u ¼ 0), recall equation (13a), that is _ For a system endowed with inputs (u  0), S v ¼ x.  v is defined in the _ Matrix H we have that Sv x. Lemma 2. 

Since S v 0, Lemma 2 guarantees Lð v; yc Þ to be quasi-definite positive, and so v; y c Þ 0 > 0 Lð v; yc Þ þ Lð

ð42Þ

_ 0pc ;1 Þ0 þ diagðB xÞ _ _ ¼ Lð~ _ 0pc ;1 Þ þ Lð~ vðxÞ; D LðxÞ vðxÞ; inequality (41) can be rewritten as: _ >0 v; yc Þ0 þ D LðxÞ Lð v; yc Þ þ Lð Provided (42), a sufficient condition for positive _ is that, for v; yc Þ0 þ  LðxÞ definiteness of Lð v; yc Þ þ Lð any vector :

A.4. Main Result

þ A01 D3 D1 ðcÞ A2 > 0

ð41Þ

By defining

Thank to assumption (19a), the system (38) implies that  ¼ 0. This, together with (35), proves that Jð; v; yc Þ > 0 for all  6¼ 0. This concludes the proof. &

_ þ A02 D3 D1 ðcÞ A1 diagðcÞ

v; yc Þ0 þ Lð~ v; 0pc ;1 Þ þ Lð~ v; 0pc ;1 Þ0 Lð v; yc Þ þ Lð _ >0 þ diagðB xÞ

       x_ x_ S 0 v þ 1 : 0pv 1 yv yv Hv

_ < 0 ðLð v; yc Þ þ Lð v; yc Þ0 Þ j0 D LðxÞj

ð43Þ

This, in turns, can be reformulated as a condition on _ For instance, sufficient the 2-norm of matrix  LðxÞ. condition such that (43) is satisfied is that:   _ 2 < min Lð v; yc þ Lð v; y c Þ 0 ð44Þ jjD LðxÞjj _ is linear in x, _ and Lð Remark that matrix  LðxÞ v; yc Þ is dependent on the values of yc and yv . Thus, there exists a strictly positive function ðyc ; yv Þ, which _ around the determines a region (for the values of x) origin where (44) holds, i.e.: _ 2 < ðyc ; yv Þ kxk

ð45Þ

Under this condition, the inequality (26) is satisfied with QðcÞ ¼ diagðcÞ. This implies, from Lemma 1, local convergence of the observer. &

We define:  v~ ¼

S Hv

~

v ¼ v  v~

x_ 0pv 1



Appendix B ð40aÞ

B. The Case of Zero Deficiency Networks ð40bÞ

where M~ denotes the pseudo-inverse matrix of M. Remark that, by construction, S v 0, and Hv v ¼ yv . _ From (40) and (29), the inequality Recall that c_ ¼ B x. (39) can be rewritten as:

In this section we deal with zero deficiency networks. Thanks to introductory concepts, lemmas and technical notions introduced in Section B.1 and in Section B.2, we prove an original statement in Section B.3. We will take advantage of this result in the proofs of Corollaries 1 and 2, given in Section B.4.

592

M. Farina and S. Bittanti

B.1. Complex Balancing and Zero Deficiency The following lemmas are related to the property of complex balancing (defined in Section 2.2) for zero deficiency networks, see [22]: Lemma 3: Given a weakly reversible zero deficiency mass-action system, let Assumption 1 hold. Then there exists an equilibrium point of equation (13b) in the positive reaction simplex S c0 , where the system is complex balanced, that is x_ ¼ 0. & Lemma 4: If a mass-action system is complex balanced at a concentration a, then it is complex balanced at all equilibrium concentrations. & B.2. Linkage Classes and Linkage Spaces In Section 2.3 we introduced the notions of linkage class and weak reversibility. Recall that l is the number of linkage classes and that vectors wi 2 R are canonical basis vectors of the complex space. Each of these vectors represents a complex of the network. One defines the linkage space Wk as:

Let Dk be the set of reactions occurring between the complexes of a linkage class Lk , i.e. Dk ¼ fwi  wj =i  jg. Let D ¼ [lk¼1 Dk . We define the set L as ð46Þ

Finally, recall the definition of the stoichiometric _ ¼ spanfn1 ;    ; nr g 2 Rn , we have space D ¼ spanfcg that: D ¼ B spanðDÞ ¼ B L

Proposition 1: Given a zero deficiency and weakly reversible network, let matrices N and S be defined as in (3a) and (7a), respectively. Then rankðSÞ ¼ rank ðNÞ ¼ s. Proof: If weak reversibility holds true, rankðSÞ ¼ dim ðDÞ ¼   l, being D defined in (46). If  ¼ 0, then also rankðNÞ ¼ s ¼   l. This proves the proposition. & Being, from Proposition 1, rankðSÞ ¼ s, ðr  sÞ linear independent vectors gi 2 R (referred to as conserved moiety vectors) such that gij s (the  entries of gi ) are integers (positive or zero) exist such that: gi S ¼ 0 ; 8 i This implies that, if r > s, by equation (13b):

Wk ¼ spanfwi =i 2 Lk g ; k ¼ 1;    ; l

L ¼ spanðDÞ:

structures called moiety-conserved cycles (see, for example [28, 33]). Such cycles take place when groups of chemicals (moieties) migrate through the network without suffering any net degradation or synthesis. We ‘‘extend’’ the notion of conserved moieties to matrix S. First of all the proof of following statement is given:

ð47Þ

B.3. Complex Balancing and Rate of Growth of x In this section we prove that there is a constant linear relation between c_ and x_ for weakly reversible zero deficiency networks. The following lemma is proved: Lemma 5: If the network has deficiency zero ( ¼ 0) and it is weakly reversible, then x_ is linearly dependent _ on c. & In order to prove Lemma 5 in Section B.3 we need intermediate results, given in the sequel. Conserved moieties for complexes. The stoichiometric matrix N shows interesting mathematical properties as the metabolic reaction network considered contains

gi x_ ¼ 0 ; 8 i ¼ 1; :::; r  s So we can identify a set of s complexes (indicated as xind ), such that the remaining r  s complex concentrations can be computed from the independent set as linear combinations. In particular, there exists a matrix S 2 R  s such that: x_ ¼ S x_ ind

ð48Þ

From equations (9) and (48): B S x_ ind ¼ c_

ð49Þ

The following proposition can be proved: Proposition 2: If the network has  ¼ 0 and it is weakly reversible, then the matrix B S has rank equal to s. Proof: As a result of Lemma 3 and Lemma 4, for zerodeficiency and weak reversible networks, c_ ¼ 0 implies x_ ind ¼ 0. This implies, from (49), that matrix B S has full column rank, i.e. rankðB S Þ ¼ s. & Proof of Lemma 5. Proof. (Lemma 5): Proposition 2 implies that there exists a pseudoinverse matrix of B DS . Let ðB DS Þ~ indicate such pseudoinverse matrix. We can write: x_ ind ¼ ðB DS Þ~ c_

ð50Þ

593

An Observer for Mass-action Reaction Networks

From (48) and (50) it follows that: x_ ¼ DS ðB DS Þ~ c_

ð51Þ

This concludes the proof.

&

B.4. Corollaries Proof (Corollary 1): Theorem 2 guarantees that observer (17) converges locally for sufficiently low rates of growth of variable xðtÞ. This is guaranteed by condition (45). From Lemma 5, under zero deficiency and weak reversibility conditions, there exists a linear _ So, for all x , there exists relation (51) between c_ and x. _ 2 < c implies that c such that the requirement kck _ 2 < x . If Assumptions (19) hold, from Theorem 2, kxk the observer (17) converges locally. & Proof (Corollary 2): If the no boundary Assumption 3 is satisfied, in addition to zero deficiency and weak reversibility, from Theorem 1 there exists a steady state c in the positive state space, globally asymptotically stable of the system if u ¼ 0. If c_ ¼ 0 then x_ ¼ 0 and condition (45) is satisfied. If Assumptions (19) hold, from Theorem 2, the observer (17) converges locally. &

Appendix C C. The case of degradation and formation reactions Proof (Corollary 3): Let B0 , S0 and S0in the system matrices in case degradation and formation are neglected. In this case, we compute Hc and Hv in order to satisfy assumptions (19). Let us consider formation and degradation separately: – If formation reactions are considered, as discussed in Section 2.4, one can consider known species synthesis rates as additional inputs. Hence, system matrices conserve their structure as in the case when formation were neglected. Hence, in order to satisfy conditions (19), Hc and Hv are the same as the case where formation is neglected. – To prove this statement we resort to the proof of Theorem 2 (see Appendix A). Note that, the main difference with the case dealt with in Theorem 2 is that the overall network is now non weakly reversible, since degradation reactions are included. So, the results of the proof not requiring weak eversibility remain unchanged. Since the only result

which requires weak reversibility is Lemma 2, we now show that such results still holds if degradation terms appear in the network topology. By including degradation, the system matrices change. We denote with B, S and Sin , these matrices in case degradation is considered. As discussed in Section 2.4, degradation can be seen as a reaction having a ‘‘virtual’’ product complex. Hence, in the network, the total number of complexes (resp. reactions) is tot ¼  þ d (resp. rtot ¼ r þ rd ), where  is the number of real complexes and d is the number of additive complexes involved in degradation. Note that the number of linkage classes increases of a number ld ¼ d  rd , while the number of ‘‘virtual’’ product complexes added is rd and the number of ‘‘real’’ degrading complexes added is ld . Without going too much in details in the structure of the matrices below, we have: 2 3 2 0 3 ð1Þ Sin S0in Sout 0rd 6 7 Sin ¼ 4 0l r Sð2Þ 5; Sout ¼ 4 0ld r 0ld rd 5 d in 0rd r Ird 0rd r 0rd rd B ¼ ½ B0

B

0nrd ;

Htot v ¼ ½ Hv

0pv rd 

while the reaction rates vector v 2 R is splitted in two components: the ‘‘classical’’ reaction rates vector v0 2 Rr , and the degradation rates vector vd 2 Rrd . 0 ð1Þ Namely, v ¼ ½v0 v0d 0 . Here S0in ; S0in 2 Rr , Sin 2 ð2Þ rd ld rd n nld 0  ; Sin 2 R ,B 2R ,B 2R . R Note that the assumption of Lemma 2 (that Sv ¼ 0) can be met only if vd ¼ 0, while v0 6¼ 0 (elements-wise) is compatible with Sv ¼ 0, since weak reversibility (for the network obtained neglecting degradation) holds. We now resort to the proof of Lemma 2 (see Section A.3). In case degradation is considered, the function Jc ðÞ is defined as in (30b) and, since vd ¼ 0, the function Jv ðÞ in (30a) can now be written as:    0 B0 0 diagðS0 0 B0 0 Þ v0 :  0 Þ0 H Jv ðÞ ¼ 0 S0in  S0out ðH v v in rþrd

Following the same arguments of Section A.3, we obtain that Jc ðÞ  0 and Jv ðÞ  0, and that the following assumptions

 0 

 0  B S0 ¼ n ; rank ¼ rank 0 Hc H0v S0in guarantee that Jc ðÞ ¼ 0 and Jv ðÞ ¼ 0 only if  ¼ 0. The remainder of the proof can be carried out following the lead of the proof of Theorem 2. &