Complexity
Volume 2018 (2018), Article ID 6879013, 20 pages
https://doi.org/10.1155/2018/6879013

Research Article

Distributed Finite-Time State Estimation of Interconnected Complex Metabolic Networks

Alfonso Sepulveda-Galvez,¹ Jesus A. Badillo-Corona,² and Isaac Chairez²

¹SEPI-UPIBI, Instituto Politécnico Nacional, 07340 Mexico City, Mexico
²Department of Bioprocesses, UPIBI, Instituto Politécnico Nacional, 07340 Mexico City, Mexico

Correspondence should be addressed to Isaac Chairez; ichairezo@gmail.com

Received 16 November 2017; Accepted 5 March 2018; Published 9 May 2018

Academic Editor: Danilo Comminiello

Copyright © 2018 Alfonso Sepulveda-Galvez et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

A set of distributed robust finite-time state observers was developed and tested to estimate the main biochemical substances in interconnected metabolic networks with complex structure. The finite-time estimator was designed by composing several supertwisting based step-by-step state observers. This segmented structure was proposed accordingly to the partition of metabolic network obtained as a result of applying the observability analysis of the model used to represent metabolic networks. The observer was developed under the assumption that a sufficient and small number of intracellular compounds can be obtained by some feasible analytic techniques. These techniques are enlisted to demonstrate the feasibility of designing the proposed observer. A set of numerical simulations was proposed to test the observer design over the hydrogen producing metabolic behavior of Escherichia coli. The numerical evaluations showed the superior performance of the observer (on recovering immeasurable state values) over classical approaches (high gain). The variations of internal metabolites inserted in the hydrogen productive metabolic networks were collected from databases. This information supplied to the observer served to validate its ability to recover the time evolution of nonmeasurable metabolites.

1. Introduction

Metabolic engineering (ME) represents one of the most relevant disciplines in bioengineering [1]. The set of methods and techniques integrated in this novel discipline seeks to optimize the metabolic circuits and regulatory mechanisms in different cells that increase their production of relevant metabolites [2]. ME is closely connected to genetic engineering and molecular biology [3] that can use complex interconnected network models with diverse structures.

ME takes into consideration the tools of applied mathematics to obtain models of the metabolic networks under study [4, 5]. The correct application of these models could save lots of resources and reduce the time to introduce productive modified cells with their remarkable secondary metabolites to the market.

There are different options to generate the model including Boolean structures, algebraic relationships, or time dependent descriptions based on ordinary or partial differential equations (ODE/PDE) [6]. This last option seems to be the most complete because the transient behavior of the cell can be captured in the ODE model. However, two natural problems arise when ODEs are used: (1) the data density needed to characterize the model is bigger than all other cases and (2) the number of parameters included in the model could be so large that existing parametric identification methods cannot be powerful enough to get a complete validation of the model [7].

The identification of uncertain parameters from time-series measurements of interesting biochemical compounds is a major aspect in system biology and ME. The majority of existing methods aimed at obtaining the parameters require the information of all compounds continuously or at least with some periodicity. On the other hand, best-fit parameter estimates are ill-posed due to issues related to data informativeness, problem formulation, and parameter sensitiveness [8, 9]. Even when the so-called canonical power-law formalism (using relative simple structures in the right-hand side of the ODE representing the metabolic network) is considered to obtain the model, the parametric identification solution remains as a complicated problem. This characteristic is emphasized when regulatory interactions among metabolites are considered. Despite the benefits of power-law formalism, the number of parameters increases if the metabolic network is described more precisely (including more metabolites and regulation interactions).

One major additional issue that must be solved to get accurate parameter values in the metabolic network representation is the necessity of measuring all the metabolites in the network. This is one of the most challenging aspects when ODEs are used to generate the model. Several experimental options have been proposed as the so-called metagenomics or using real-time polymerase chain reaction (PCR). Nevertheless, the expensiveness of applying these techniques may limit their application on characterizing metabolic networks. The number of experiments needed to characterize metabolic networks accurately is usually large. Nevertheless, in actual metabolic networks, only a small fraction of intracellular metabolites can be directly measured, and therefore initial conditions should be also estimated. Undoubtedly, in metabolic networks, lack of experimental data is unavoidable. This condition compromises the accurateness of parameters estimated by any feasible method [10].

State estimation in metabolic engineering has become an option to recover the time variation of compounds in metabolic networks. Different options have been proposed to solve the reconstruction of the metabolites concentration over time. However, the complexity of these networks limits the application of global observers for the entire network. One popular strategy recently explored is the divide-and-conquer scheme where a set of low-order state observes are running in parallel. Each of these observers is applied on a certain section of the network where observability property holds [11–13].

State observers to recover immeasurable information from biological and biotechnological systems have been developed for many years [14]. The application of these observers can provide information that can be eventually used to regulate the metabolic network and, in consequence, optimize the production of some metabolites. In this sense, a state observer is known as a software sensor. That is, a suitable and well-designed state estimator can be used as an accurate artificial sensor of some specific variables. The artificial measurements provided by the observer correspond to variables that cannot be measured online or their measuring cost is relatively high [15].

This study describes the design of a so-called step-by-step decentralized observer to estimate the unknown states of the selected metabolic network. This state estimator is composed of a set of robust high-order sliding mode differentiators. In particular, the decentralized observer is applied on a system representing hydrogen () production by a strain of Escherichia coli.

Notice that, in this study, no general characteristics of metabolic networks are used for the construction of the observers or for proposing a method to estimate their state variables. We only considered a simplified metabolic network that cannot generalize all the constraints, restrictions, and characteristics of general metabolic networks. Any possible solution for such problem requires a deeper understanding of internal interactions between metabolites, enzymes, and genes. However, a possible solution of that problem is beyond the scope of this article.

2. Notation

This section introduces the nomenclature used in this manuscript. Acronyms of all substances considered to construct the model of metabolic network are listed as follows.

 : Biomass Glc: Glucose G6p: Glucose 6-phosphate 2gp: 2-Phospho-D-glycerate Pep: Phosphoenolpyruvate Oxa: Oxalacetate Mal: Malate Suc: Succinate Cit: Citrate -kg: -Ketoglutarate Pyr: Pyruvate Acoa: Acetyl coenzyme A Aad: Acetaldehyde EtOH: Ethanol Acp: Acetylphosphate Acet: Acetate For: Formate Lac: Lactate H₂: Hydrogen CO₂: Carbon dioxide.

3. Hydrogen Producing Metabolic Network of E. coli

The model used in this study was built using 18 ODEs that represent the dynamic behavior of all metabolites involved in the metabolic pathway presented in Figure 1. Notice that this system was selected just as an example that demonstrates the application of the state estimation technique proposed in this study. The first ODE of the model represents the production of biomass and the rate of the reaction is modeled as Monod kinetics containing two different parameters: and . The subsequent reactions correspond to intracellular and extracellular metabolites of the metabolic network. These equations consider the production and consumption of each metabolite due to the action of the enzymes that catalyze the reactions. Each enzyme has a specific reaction rate depending on the concentration of its substrate, product, or any other metabolite that participates in the regulation mechanisms and can take a positive sign to express production or a negative sign to express consumption. The cell takes up Glc from the culture media by the action of the PTS system. This system is activated by the presence of Pep and takes Glc as substrate forming G6p as a product. According to [17], an activated mechanism can be modeled as a bisubstrate reaction (a second-order one).

Figure 1: Metabolic pathway of hydrogen production in Escherichia coli.

The conversion of G6p to Fbp is made by the enzyme pfkA. This reaction is inhibited by the presence of Pep and it is modeled as a bisubstrate reaction due to the presence of ADP and taking into account Pep as an inhibitor.

Enzymes gpmA and eno catalyze the transformation of Fbp to 2pg and 2pg to Pep, respectively. These reactions are modeled as a Michaelian first-order system.

The variable Pep is used in two different reactions, one of them to produce Oxa by the enzyme ppc that is repressed by the presence of Mal and activated by Acoa and Fbp. This reaction rate is modeled as a trisubstrate reaction [17] due to the action of Acoa and Fbp as activators and taking into account Mal concentration as an inhibitor. Pep is also converted to Pyr by the enzyme pyk that is activated by Fbp whereby it is modeled as a bisubstrate reaction.

The term Oxa is converted to two different products by two different reactions. In this reaction, the involved enzyme is mdh. This enzyme uses Oxa as a substrate and the products of this reaction are Mal and NADH. The second reaction is catalyzed by the enzyme gltA. This reaction takes Oxa and Acoa as substrates to produce Cit and it is inhibited by the concentration of -kg.

The variable Mal is converted to Suc by the action of the enzyme frdA that is inhibited by Oxa concentration. Cit is converted into two products by the action of the enzyme icd. Pyr is used as a substrate in two different reactions; the first one is catalyzed by the enzyme ldhA to produce Lac that is in turn excreted to the media culture. In the second reaction, Pyr is converted to Acoa and For by the enzyme pflB. Formate is used by the complex Fhc to obtain H₂ and CO₂. Also, the enzyme focA excretes For to the culture media. Acoa is also used by the enzymes pta and adhE. Enzyme pta is activated by Pep to produce Acp that is converted to Acet by the enzyme ackA and released to the media culture. Enzyme adhE converts Acoa to Acd that is eventually converted into EtOH.

The proposed modeling strategy yields the mathematical model of the metabolic pathway focused on the H₂ production by Escherichia coli. The model is presented in (1). This model contains all the equations considered as a solution of the thermodynamic study with the set of parameters gotten from bibliography sources. The model proposed in this study has not been presented before and it was constructed using the arguments described above.

The set of parameters used in the model presented in (1) is detailed in Table 1. These parameters were obtained from bibliographic references as shown at the bottom of the table. Many of them have been experimentally confirmed in different experimental analytic forms. The initial conditions used in the numerical simulation results were also taken from a similar set of published results.

Table 1: Parameter values used to simulate the metabolic network.

Parameter
Real value (h⁻¹)
Parameter
Real value (mM)

Parameter
Real value (mMh⁻¹)
Parameter
Real value (mM)

Taken from [16]. Taken from EcoCyc database.

4. Segmentation of Metabolic Network

The entire metabolic network was divided into 7 different subsystems. The rules to divide the network were proposed in order to get simpler models that satisfy the observability condition proposed in [18]. Then, a set of seven different output-input pairs were selected in order to satisfy the aforementioned conditions in each subsystem. These conditions were obtained from the structural analysis of networks, including the formation of cycles and bifurcations. Based on these rules, Table 2 contains the information of which states were selected as input and output in each subsystem.

Table 2: Distribution of input-output pairs in the different subsystems of the metabolic network.

Subsystem
Input variable
Output variable

1
Glc
Pep

2
Pep
Suc

3
AcoA
-kg

4
AcoA
EtOH

5
AcoA
Acet

6
Pyr
H₂

7
Pyr
Lac

Each of the input-output pairs was selected in order to use the canonical Brunovskii transformation. The selection of the corresponding output variable provides the full relative degree for each subsystem according to the procedure described in [19]. Therefore, using the results given in the same reference, there exists a nonlinear transformation , where

where is the order of each subsystem . By assumption, all the time derivatives are different from zero for all and bounded at least to the th derivative. If the metabolic network has no regulation feedback or it can be neglected because of its relative impact on kinetic equations of each metabolite, then the observer proposed in this study can be applied. In consequence, the new variables obey the following dynamics:

where the nonlinear functions and can be calculated as a solution of a regular differentiation method of the corresponding output variable. The functions are obtained by the presence of feedback regulations. However, all of them are measurable because they depend on those metabolites selected as inputs or outputs of the submodels described in Table 3. The variables marked with represent all the states of subsystems different from the th one.

Table 3: Real and estimated parameters and relative error between them.

Parameter
Real value (h⁻¹)
Estimated value (h⁻¹)
Relative error (%)
Parameter
Real value (mM)
Estimated value (mM)
Relative error (%)

1.002

1.666
0.96

1.0

Parameter
Real value (mMh⁻¹)
Estimated value (mMh⁻¹)
Relative error (%)
Parameter
Real value (mM)
Estimated value (mM)
Relative error (%)

5.0

5.0

4.98

2.92

4.97

5.0

5.0

0.36

5.0

5.0

0

5.0

0

0

0

0

5

5.0

5.0

5

5

5

4.97

4.98

4.97

1.4

5.0

1.42

4.77

2.31

5.0

1.84

4.97

5.0

4.63

2.2

*All the corresponding real values were obtained from [16]. **The real parameters came from [6].

In all the cases considered in this study, all the functions were different from zero. Therefore, the relative degree is completed and well defined in all the selected subsystems. The possible situations where the relative degree may be less than are never attainable in the system proposed to represent the metabolic network. Based on this condition, each subsystem is observable with the corresponding input-output pairs selected as in Table 3. Notice that full relative degree is not a necessary condition because only the observability restriction is needed within each subsystem [20]. Therefore, a simpler version of the nonlinear transformation can be proposed but this is a matter of further research.

4.1. Example of Transformation Procedure

Let us consider the first subsystem characterized using Glc as input and Pep as output. So, the nonlinear transformation between the set of states Pep, 2pg, Fbp, G6p, and Glc and the corresponding is defined by

where is Glc.

4.2. Basic Assumptions

The corresponding calculus of nonlinear functions and yields the following assumptions.

Assumption 1. The function associated with the control action is bounded and does not change sign:

where and   are two positive constants. The input signal must fulfill the following.

Assumption 2. The control action belongs to the following admissible set:

Assumption 3. All the functions satisfy the following inequality:

where are positive scalars.

Assumption 4. By the nature of elements included in the nonlinear transformation , all the variables are bounded as follows:

The upper bounds described in the previous assumption can be obtained by the time evolution of metabolites or they can be theoretically calculated by the mass balance executed on the specific metabolic network considered in this study as presented in [10].

5. Finite-Time Convergent Observer of Segmented Metabolic Network

The set of observers proposed in this study satisfies the following differential equation:

with the estimation error satisfying the following equation:

The set of variables represents the corresponding estimated trajectories of . The variables serve as auxiliary (intermediate) estimation of velocity for each variable . In particular, . The observer gains represented by , and are scalars that must be adjusted to force the convergence of the observer trajectories to the states of the uncertain system.

The gains and belong to the set

The indicator function fulfills the following definition:

The switching time is found as a result of the fixed-time converge obtained for the observer in the following section. The nonlinear functions and were designed in agreement with the proposal given in [21]. Then, the following structures were considered in each observer design:

Parameters represent extra gains that should be adjusted according to the results presented in the main theorem of this work. For the purposes of this study, the function was used as follows:

One must recall that the solution of (9) is understood in the sense of Filippov [22]. The observer proposed in this study provides the so-called fixed-time stability for the origin of the estimation error. To clarify this concept, we introduced the following definition [23].

Definition 5. The origin of

where and , is said to be (a) globally finite-time stable if it is globally asymptotically stable and any solution reaches the equilibrium point at some finite time () or (b) globally fixed-time stable if it is globally finite-time stable and the convergence time is bounded, that is, fulfilling

Based on the previous definition, the convergence for each individual observer has been already proven in [24]. Nevertheless, the solution of the observer proposed in this study requires a new theoretical result described in the following proposition.

Proposition 6. Suppose that assumptions (5) and (6) are valid; then, the estimation error (10) is uniformly exact convergent to the origin with a fixed-time given by

if there exist some positive scalars such that the following pair of LMIs have at least one positive definite solution of each ,  , where each matrix is the positive definite solution of

Proof. Proving the convergence (in uniform exact sense) of the states of each observer given in (9) to the variables of (3) requires the application of two Lyapunov-like candidate functions [21]. The first function is used to demonstrate the global and exact convergence of (9) to the trajectories of (3), but not the uniformity with respect to the initial conditions. The second function is proposed to show the robust uniform convergence of the estimation error trajectories.
The convergence of each state observation error is obtained in steps.
Step 1. Assume that and . Let us introduce a set of auxiliary variables defined as (). The error between the states of each (9) and the corresponding transformed system (3) is characterized by and Their dynamics within the period of time obey

The first couple of equations in (18) coincide with the so-called modified supertwisting algorithm. The rest of the dynamics do not grow because are all identically zero, considering the nature of the functions and
To prove the fixed time of , two Lyapunov functions are needed [21]: the first one provides the finite-time stability of but it cannot be used to prove the fixed-time stability. The second Lyapunov function is used to ensure the fixed-time stability of but considering that is not so far from the origin. The first energetic function used in is defined as :

The second energetic function used in the same period of time is defined as

where and are two positive scalars. According to the result presented in [21], the combination of both energetic functions can be used for proving that ,  , under the assumption of , with . Therefore, if , the following identities hold:

Step 2. Within the period of time , the following identity is true: ; then, the dynamical evolution of and is given by

Once again, the dynamic equations of and coincide with the so-called modified supertwisting algorithm. The rest of the dynamics do not grow because are identically zero. A new combination of two energetic functions named and defined by

are used in this part of the proof. Once again, these two functions can be used for proving that ,  , and therefore

Step . The last step of the observer design included the accurate design of to stabilize

The last part of the proof is based on a Lyapunov-candidate function defined as follows:

The regular analysis based on Lyapunov technique yields

The substitution of in the previous differential equation leads to

The application of the Cauchy-Schwartz inequality provides the following upper bound for :

By virtue of the continuity of all functions and the bounds considered in Assumption 2, one gets

where . Hence, and Then, considering the value of presented in (5), one gets

Following a regular procedure [25], one can show that

6. Parametric Characterization of Metabolic Network

It is evident that each differential equation of the model presented above satisfies

where is any variable in the model presented in (1). The entire set of variables is lumped in the supertwisting algorithm vector and is the number of variables. Then, represents all the variables involved in the metabolic network under study. Also, it is assumed that all the functions and are linearly parametrized by and , respectively. Nonlinear functions ,   satisfy

Therefore, a couple of nonlinear functions and exist, such that

The nonlinear functions and are introduced to define the so-called residual representation of and . In the unrealistic case when the derivative of the whole vector can be measured, the ODE described in (33) can be alternatively represented as

This form obeys the regular linear regression form. The solution for the parametric identification can be obtained as in [26]:

Table 3 demonstrates the real and estimated values of all the parameters calculated by the method proposed in this study.

7. Implementation Issues

7.1. Effect of Sampling the Metabolites Concentrations

The set of finite-time state estimators proposed in this study requires a special implementation scheme. This section details how the set of distributed finite-time state estimators was developed.

The keystone technical fact to implement the observer proposed in this study is the possibility of measuring all the metabolites considered as input or output in the segmentation process. Nowadays, there are several methods to obtain the actual concentrations of key internal metabolites in the cell.

For example, the study presented in [24] reports a detailed list of experimental methods to detect all the metabolites considered either as input or as output of each subsystem proposed in this study. Despite the fact that the observer suggested in this study demands the continuous monitoring of metabolites concentrations, it has been proven that the sampling process decreases the quality of estimation. Nevertheless, this reduction is not severe if the sampling period is small enough. For example, in the case of Escherichia coli cultures, this time can be as small as 2 minutes. Then, the observer proposed in this study still works on the more realistic case of sampling the information of metabolites concentrations.

The estimation method for substances considered in the metabolite network assumes that their concentrations can be obtained continuously which may seem as an unrealistic assumption. Nevertheless, the effect of discrete sampling in the output/input information injected in each sliding mode observer has been previously analyzed in different studies such as [27–29]. All of these studies explained the losses in the estimation quality as a consequence of the sampling process. According to [30], each estimation step performed by the observer proposed in this study introduces an estimation error which is proportional to the square of sampling time and proportional to the upper bound of uncertainties affecting the supertwisting dynamics. Notice that this sampling induced deviation prevents the finite-time convergence to the origin of the estimation error at each step. In consequence, the accumulated error after -estimation steps cannot be estimated directly as a proportional value of the estimation error obtained after the first step of differentiation. Indeed, this is an open problem in the implementation of robust step-by-step differentiators under sampling process. Therefore, if metabolites are measured discretely, the effect is a losing of estimation quality with a bound defined in previous studies. This is a natural consequence that has been observed not only in the case of sliding mode observers but also in high-gain state estimators and some others.

7.2. Effect of Noisy Measures of Metabolites Concentrations

A second relevant assumption regarding the observer developed in this study is the noise presence on the measured signal. Noise plays a major role in the estimation quality obtained by any sliding mode based observer.

The set of observers proposed in this study assumed that is measured free of noises. Under this assumption, the estimate of the state is achieved in fixed time. However, if the output signal is measured affected by noise , that is, with being the noise signal, the result presented in the theorem and using the results presented in [31] showed that the estimation error is input to state stable if the noise is considered as the external unknown input. On the other hand, STA is particularly sensitive to the noise effect. This characteristic can be attenuated by a correct selection of gains. The study presented in [32] reviewed deeply this issue. The first theorem in that study claims that the ultimate bound for denoted by is the largest solution of

Based on the value obtained in (38), the ultimate bound of , namely, , is given by

where the notation denotes the th element of a two-dimensional matrix. This study did not consider the presence of noises. However, if experimental information is injected to the observer structure, then the method described in (38) can be applied.

8. Numerical Simulation

This section describes the solution obtained when the proposed observer was evaluated numerically. The observers were simulated using the trajectories of model (1) with the parameters collected from different published results regarding kinetic characterization of Escherichia coli. These parameters were obtained in different studies. A summary of the parameters is presented in Table 2 where the references consulted for simulation purposes are detailed.

The production of was selected as the target metabolic network considering that research on renewable fuels has gained importance due to the decline in international reserves of fossil fuels as well as their climatic and economic effects [33, 34]. Among the investigated renewable fuels, has some advantages; for example, it has high energy efficiency (122 KJg⁻¹) and its application produces only water as a byproduct and can be obtained from cheap carbon sources [35–37].

production can be obtained as a result of microorganisms biological activity. Among other advantages of these methods, the possibility of using different nonconventional sources (even waste resources) as a carbon source for microorganisms seems to be a relevant factor when technical-economical balance is a key condition. Moreover, many of these processes may occur under normal conditions of pressure and temperature [35, 38].

Biological hydrogen production can be made with different microorganisms. One of the most popular is Escherichia coli, a Gram-negative facultative anaerobic bacterium, which is able to produce by mixed acid fermentation with a theoretical maximum production of 2 moles of by 1 mole of glucose [39]. Main efforts to improve production have focused on (a) identifying and characterizing the effect of oxygen-tolerant hydrogenases, (b) improving molar yields, (c) developing efficient separation techniques from bioreactor head-space, and (d) modifying metabolic pathways response by genetic engineering [40].

The optimization of production can be obtained by regulating substrates used to feed Escherichia coli. These substrates may not serve only as carbon or nitrogen sources but as genetic regulator. The expensiveness of this evaluation can compromise the economical yield produced by the production of Escherichia coli with respect to the resources invested in final optimization. The accurate model of the H₂ production by E. coli can be used to test diverse genetic or metabolic strategies. This testing scheme can be used to reduce the final production costs. This is a main motivation of this study.

The parameter values were taken as they were presented in all the studies described in Table 4. The numerical evaluation of transformation for system (1) presented in this study was simulated in Matlab/Simulink. The simulation was executed using a fixed step numerical integration algorithm (ODE8-Dormand Price method) with an integration step of 0.0001 hours. This value has been recommended when simulations include some variants of sliding mode systems [31]. In order to offer a fair comparison between STA and high-gain observers, a similar gain tuning process was considered including similar time processing. The gains in the STA solution only use the values predicted by the set gain included in the main theorem presented in this study. The parameters included in the high-gain observer satisfied the tuning rules presented in [41]. The set of initial conditions used in simulation appears in Table 4.

Table 4: Initial conditions on space and original space.

Subsystem 1

Pep

pg

G6p

Subsystem 2

Suc

Fum

Mal

Subsystem 3

kg

Iso

Cit

Subsystem 5

EtOH

Aad

Subsystem 6

Acet

Acp

Subsystem 7

H₂

For

Figure 2 shows the performance of the first observer in the -coordinates for the first subsystem. The trajectories depicted in this figure were obtained as a solution of the set of discontinuous systems described in (3). These figures show also the step-by-step convergence of all the states. The first state (Figure 2(a)) converged before (Figure 2(b)) and it did it before (Figure 2(c)). The convergence of all states is achieved in finite time and the estimation error converged before 0.5 h. This value coincides with the predicted time proposed in the theorem and it is considered to be relatively low with respect to the period of simulation (30.0 h). The same set of figures shows the state estimation attained as a solution of the first high-gain observer applied on the first subsystem. At first sight, no evident differences are detected between the high-gain observer and STA trajectories. Nevertheless, the smaller differences observed within the first two hours of reaction are enough to induce high frequency oscillations in the states of the high-gain observer, especially after applying the inverse transformation.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: Simultaneous visualization for the states of the transformed coordinates in the first subsystem, the corresponding high-gain observer, and the STA-based estimator: (a) , (b) , and (c) . (d)–(f) depict the variation of the concentrations for the compounds (d) , (e) , and (f) as well as their estimates produced by the inverse transformation applied to both the high-gain and the STA-based observer.

The convergence of the estimated states in the -coordinates forces the corresponding fixed-time estimation of all the states in the first subsystem that includes the effect of substrate. These values are obtained as a solution of the inverse transformation presented above. Figures 2(d), 2(e), and 2(f) compare the time evolution of G6p, pg, and Pep. Notice the presence of a high amplitude over estimation in the estimation of G6p.

The graphical information included in Figure 2 may be insufficient to justify the superior performance of STA-based observer. Figure 3 depicts the variation of the same states shown in Figure 2. However, the time period of state evolution was constrained to the first 5 hours of reaction. These closer views to the earlier stages of the reaction show how the finite-time convergence is attained when only two key compounds are measured (input-output). Indeed, the number of estimated states where high frequency oscillations and high amplitude transient states were detected was smaller when the STA observer was applied to recover the nonmeasurable states of the first subsystem.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Detailed simultaneous visualization for the states of the transformed coordinates in the first subsystem, the corresponding high-gain observer, and the STA-based estimator: (a) , (b) , and (c) during the first 5 hours. (d)–(f) depict the detailed variation of the concentrations for the compounds (d) , (e) , and (f) as well as their estimates produced by the inverse transformation applied to both the high-gain and the STA-based observer during the first 5 hours.

Figure 4 depicts the time evolution of the trajectories included in the third subsystem. The transformed states ,  , and corresponding to citrate, isocitrate, and -kg as well as their estimates ,  , and in the transformed states are shown in Figures 4(a), 4(b), and 4(c). This set of figures evidences the remarkable transient behavior of the observer proposed in this study. No oscillations of high frequency and amplitude were obtained when the STA-based observer was considered. Nevertheless, if a high-gain observer is considered, unacceptable oscillations of the estimated states are gotten for the estimated states. These oscillations were obtained even after applying different strategies to adjust the observer gains. Actually, the state estimation performance is deficient in the early stages of the reaction. This can be considered relevant because the duplication period of Escherichia coli is around 45 minutes. The corresponding states obtained after applying the corresponding inverse transformation and the estimated states of Cit,  IsoCit, and -kg matched the actual trajectories of the corresponding states from model 1. One may notice that there is no evident difference between the three pairs of trajectories when the STA observer was applied. Despite the measurable variable considered to adjust the observer performance, the observer proposed in this study showed a lower level of oscillations around the true trajectories of the metabolic network and larger convergence time.

(a)

(b)

(c)

(d)

Figure 4: Simultaneous visualization for the states of the transformed coordinates in the third subsystem, the corresponding high-gain observer, and the STA-based estimator: (a) and (b) . (c) and (d) depict the variation of the concentrations for the compounds (c) Aad and (d) EtOH as well as their estimates produced by the inverse transformation applied to both the high-gain and the STA-based observer.

Figure 4 justifies how the performance of STA-based observer was better than the one obtained by the high-gain state estimation. Figure 5 shows the time variation of the transformed estimated states. The variation of the same states shown in Figure 4 is detailed within the first 5 hours of reaction. Once more, these closer views to the initial period of the reaction show the finite-time convergence of the STA-based observer. Notice that states of the STA-based observer approach the actual states while the high-gain version actually touches the real trajectories but they did not remain over the trajectory as shown in 4 which shows a longer time period.

(a)

(b)

(c)

(d)

Figure 5: Detailed simultaneous visualization for the states of the transformed coordinates in the third subsystem, the corresponding high-gain observer, and the STA-based estimator: (a) and (b) during the first 5 hours. (c) and (d) depict the detailed variation of the concentrations for the compounds (c) Aad and (d) EtOH as well as their estimates produced by the inverse transformation applied to both the high-gain and the STA-based observer during the first 5 hours.

Figure 6 details the dynamical behavior of transformed trajectories considered as part of the seventh subsystem. Both transformed states and (that can be transformed into formate and hydrogen) as well as their estimates and in the transformed states are shown in Figures 6(a) and 6(b). Figures 6(a) and 6(b) contain the time evolution of and , respectively. This particular subsystem was elected to show its trajectories because it is the most simple among the seven ones considered to describe the metabolic network studied in this article.

(a)

(b)

(c)

(d)

Figure 6: Simultaneous visualization for the states of the transformed coordinates in the seventh subsystem, the corresponding high-gain observer, and the STA-based estimator: (a) and (b) . (c) and (d) depict the variation of the concentrations for the compounds (c) For and (d) H₂ as well as their estimates produced by the inverse transformation applied to both the high-gain and the STA-based observer.

The actual concentrations of H₂ and formate obtained after applying the corresponding inverse transformation are also depicted in Figures 6(c) and 6(d). The robustness of the observer proposed in this study confirms that no evident difference between both trajectories was obtained. Despite the subsystem analyzed from the metabolic network, the observer proposed in this study showed that true trajectories of the metabolic network can be reproduced by the distributed observer based on the sequential application of the high-order sliding mode.

An additional figure (Figure 7) was included to describe the variation of and as well as the time evolution of formate and hydrogen. Even when these trajectories correspond to the most simple subsystem, the superior performance of the STA-based observer is conserved. Even when both observers converge in less than one hour, the STA-based observer did it within the first half hour.

(a)

(b)

(c)

(d)

Figure 7: Detailed simultaneous visualization for the states of the transformed coordinates in the seventh subsystem, the corresponding high-gain observer, and the STA-based estimator: (a) and (b) during the first 5 hours. (c) and (d) depict the detailed variation of the concentrations for the compounds (c) For and (d) H₂ as well as their estimates produced by the inverse transformation applied to both the high-gain and the STA-based observer during the first 5 hours.

To characterize the impact of the noise and sampling process over the estimation quality for the observer developed in this study, a specific set of simulations were evaluated. The study of the noise impact was developed by aggregating a source of noise over the measurable information obtained in each subsystem. A bounded random signal generator was added to the available information with a bound of 0.02 mM. This value was proposed arbitrarily. The entire set of simulations was executed to estimate all the states in the metabolic model. Then, the upper bound for the noise was increased 10 times with an increment of 0.02 mM. The mean square value for the estimation error was calculated for the entire set of experiments. The value calculated for the mean square was compared with the upper bound proposed for the noise. Figure 8(a) presents this evaluation.

(a)

(b)

Figure 8: Evaluation of the noise (a) and sampling process (b) over the estimation quality.

The evaluation of the sampling process used a similar strategy to the one presented above (used to evaluate the effect of noise). Therefore, the sampling information was injected in the observer structure considering subsequent increments of the sampling period. The increment was of 0.2 hours (12 minutes). A set of ten increments was evaluated. Figure 8(b) depicts the corresponding evaluation of sampling effect.

In both cases, there is a loss of estimation quality (measured in terms of the mean square error increment) as a consequence of two aspects: the sampling period and the presence of noise.

The type of state estimator proposed in this study has proven to be a useful tool to reconstruct the time evolution of compounds included in metabolic networks. This fact is supported by the evident superior performance obtained when the STA-based observer was compared with a standard method of estimation such as the high-gain scheme. Even when the assumption of measuring some key compounds in the metabolic network may seem difficult to fulfill, there are a plenty of results reporting the application of diverse analytic methods that can be used to solve this part of the problem.

The finite-time convergence of the estimation error obtained in all the subsystems also provides the additional benefit of using the estimated information to solve the parametric characterization of the mathematical model. This condition has been also tested in this study where the information of the estimated states is injected into the linear regressor that is used in the least mean square method.

9. Conclusions

This study has developed a new class of state estimators to reconstruct internal metabolites in the metabolic network focused on producing hydrogen in Escherichia coli. The estimator satisfied the step-by-step structure based on the sequential application of the uniform supertwisting algorithm. The method proposed in this study used a segmentation process of the entire metabolic network. The segmentation was enforced to keep the observability condition in each subsystem. A specific step-by-step supertwisting estimator was used to estimate the metabolites in each subsystem. This distributed scheme was presented for the first time in this study. The application of this estimation method was tested in numerical simulation with better solutions than the ones obtained by high-gain observers. The simulation used a mathematical model that considered the key section of Escherichia coli’s metabolic network focused on producing its main secondary metabolites. A set of numerical simulations were executed to prove the effectiveness of the proposed STA-based observer. The comparison of the estimates obtained by the STA-based observer with respect to the regular high-gain version of the state estimator supported also the more effective performance of the suggested observer. Three of the seven subsystems were numerically reported where the estimates of the STA were closer to the states of the metabolic network. The rest of the four subsystems showed similar results but they were not included in the manuscript because of article length restrictions. The estimation proposed in this study can be easily extended to different metabolic networks in diverse microorganisms. The only condition needed to apply the technique developed in this study is to keep the observability conditions in each subsystem obtained from the segmentation process.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.

References

R. Liu, M. C. Bassalo, R. I. Zeitoun, and R. T. Gill, “Genome scale engineering techniques for metabolic engineering,” Metabolic Engineering, vol. 32, pp. 143–154, 2015.

D. Jullesson, F. David, B. Pfleger, and J. Nielsen, “Impact of synthetic biology and metabolic engineering on industrial production of fine chemicals,” Biotechnology Advances, vol. 33, no. 7, pp. 1395–1402, 2015.

A. Fernández-Castané, T. Fehér, P. Carbonell, C. Pauthenier, and J.-L. Faulon, “Computer-aided design for metabolic engineering,” Journal of Biotechnology, vol. 192, part B, pp. 302–313, 2014.

Y. Zheng, C.-D. Yang, J.-W. Yeh, and S.-S. Jang, “Metabolic regulatory network optimization using an information guided genetic algorithm approach,” Computer Aided Chemical Engineering, vol. 21, pp. 1705–1710, 2006.

Z. Sapci and J. Morken, “The effect of algae species on biodiesel and biogas production observed by using a data model combines algae cultivation with an anaerobic digestion (ACAD) and a biodiesel process,” Energy Conversion and Management, vol. 79, pp. 519–524, 2014.

T. A. A. Kadir, A. A. Mannan, A. M. Kierzek, J. McFadden, and K. Shimizu, “Modeling and simulation of the main metabolism in Escherichia coli and its several single-gene knockout mutants with experimental verification,” Microbial Cell Factories, vol. 9, article 88, 2010.

M. Ashyraliyev, Y. Fomekong-Nanfack, J. A. Kaandorp, and J. G. Blom, “Systems biology: Parameter estimation for biochemical models,” FEBS Journal, vol. 276, no. 4, pp. 886–902, 2009.

I.-C. Chou and E. O. Voit, “Recent developments in parameter estimation and structure identification of biochemical and genomic systems,” Mathematical Biosciences, vol. 219, no. 2, pp. 57–83, 2009.

M. Ederer, S. Steinsiek, S. Stagge et al., “A mathematical model of metabolism and regulation provides a systems-level view of how Escherichia coli responds to oxygen,” Frontiers in Microbiology, vol. 5, article 124, 2014.

K. Pontes Freitas Alberton, A. L. Alberton, J. A. Di Maggio, V. G. Estrada, M. S. Díaz, and A. R. Secchi, “Simultaneous parameters identifiability and estimation of an E. coli metabolic network model,” BioMed Research International, vol. 2015, Article ID 454765, pp. 1–21, 2015.

Y.-Y. Liu, J.-J. Slotine, and A.-L. Barabási, “Observability of complex systems,” Proceedings of the National Acadamy of Sciences of the United States of America, vol. 110, no. 7, pp. 2460–2465, 2013.

A. J. Whalen, S. N. Brennan, T. D. Sauer, and S. J. Schiff, “Observability and controllability of nonlinear networks: The role of symmetry,” Physical Review X, vol. 5, no. 1, Article ID 011005, 2015.

F. Cacace, A. Germani, and P. Palumbo, “The state observer as a tool for the estimation of gene expression,” Journal of Mathematical Analysis and Applications, vol. 391, no. 2, pp. 382–396, 2012.

J.-P. Gauthier, H. Hammouri, and S. Othman, “A simple observer for nonlinear systems applications to bioreactors,” Institute of Electrical and Electronics Engineers Transactions on Automatic Control, vol. 37, no. 6, pp. 875–880, 1992.

J. A. Moreno, J. Alvarez, E. Rocha-Cózatl, and J. Diaz-Salgado, “Super-twisting observer-based output feedback control of a class of continuous exothermic chemical reactors,” IFAC Proceedings Volumes, vol. 43, no. 5, pp. 727–732, 2010.

K. R. A. Castro, Determinación de los flujos metabólicos en la producción de hidrógeno [Ph.D. thesis], Unidad Profesional Interdisciplinaria de Biotecnologa-IPN, 2011.

V. Leskovac, Comprehensive Enzyme Kinetics, Springer, New York, NY, US, 2003, http://books.google.com.mx/books?id=98mDseisahEC.

W. M. Haddad and V. Chellaboina, Nonlinear Dynamical Systems and Control. A Lyapunov-Based Approach, Princeton University Press, Princeton, NJ, USA, 2008.

A. Isidori, Nonlinear Control Systems, vol. 1 of Communications and Control Engineering Series, Springer, Berlin, Germany, 1995.

J.-P. Gauthier and G. Bornard, “Observability for any of a class of nonlinear systems,” IEEE Transactions on Automatic Control, vol. 26, no. 4, pp. 922–926, 1981.

E. Cruz-Zavala, J. A. Moreno, and L. . Fridman, “Uniform robust exact differentiator,” Institute of Electrical and Electronics Engineers Transactions on Automatic Control, vol. 56, no. 11, pp. 2727–2733, 2011.

A. F. Filippov, Differential Equations with Discontinuous Right-Hand Side, Kluwer, Dordrecht, The Netherlands, 1988.

A. Polyakov, “Nonlinear feedback design for fixed-time stabilization of linear control systems,” IEEE Transactions on Automatic Control, vol. 57, no. 8, pp. 2106–2110, 2012.

N. Martínez-Fonseca, I. Chairez, and A. Poznyak, “Uniform step-by-step observer for aerobic bioreactor based on super-twisting algorithm,” Bioprocess and Biosystems Engineering, vol. 37, no. 12, pp. 2493–2503, 2014.

V. Utkin, J. Guldner, and J. Shi, “Sliding mode control in electro-mechanical systems,” in Automation and Control Engineering, 2nd edition, 2009.

L. Ljung, System Identification: Theory for the User, Prentice Hall, 2nd edition, 1999, http://www.worldcat.org/isbn/0136566952.

I. Salgado, L. Fridman, O. Camacho, and I. Chairez, “Discrete time supper-twisting observer for 2n dimensional systems,” in Proceedings of the 8th International Conference on Electrical Engineering, Computing Science and Automatic Control, CCE '11, Merida City, Mexico, October 2011.

I. Salgado, I. Chairez, B. Bandyopadhyay, L. Fridman, and O. Camacho, “Discrete-time non-linear state observer based on a super twisting-like algorithm,” IET Control Theory & Applications, vol. 8, no. 10, pp. 803–812, 2014.

N. K. Sharma and S. Janardhanan, “Discrete higher order sliding mode: concept to validation,” IET Control Theory & Applications, vol. 11, no. 8, pp. 1098–1103, 2017.

M. Livne and A. Levant, “Proper discretization of homogeneous differentiators,” Automatica, vol. 50, no. 8, pp. 2007–2014, 2014.

A. Levant, “Robust exact differentiation via sliding mode technique,” Automatica, vol. 34, no. 3, pp. 379–384, 1998.

M. T. Angulo, J. A. Moreno, and L. Fridman, “The differentiation error of noisy signals using the Generalized Super-Twisting differentiator,” in Proceedings of the 51st IEEE Conference on Decision and Control, CDC '12, pp. 7383–7388, Maui, Hawaii, USA, December 2012.

J. Mathews, Q. Li, and G. Wang, “Characterization of hydrogen production by engineered Escherichia coli strains using rich defined media,” Biotechnology and Bioprocess Engineering, vol. 15, no. 4, pp. 686–695, 2010.

P. C. Hallenbeck, “Fermentative hydrogen production: Principles, progress, and prognosis,” International Journal of Hydrogen Energy, vol. 34, no. 17, pp. 7379–7389, 2009.

T. Maeda, V. Sanchez-Torres, and T. K. Wood, “Hydrogen production by recombinant Escherichia coli strains,” Microbial Biotechnology, vol. 5, no. 2, pp. 214–225, 2012.

M. Z. Mohd Yusoff, Y. Hashiguchi, T. Maeda, and T. K. Wood, “Four products from Escherichia coli pseudogenes increase hydrogen production,” Biochemical and Biophysical Research Communications, vol. 439, no. 4, pp. 576–579, 2013.

S. K. S. Patel and V. C. Kalia, “Integrative biological hydrogen production: an overview,” Indian Journal of Microbiology, vol. 53, no. 1, pp. 3–10, 2013.

L. M. Rosales-Colunga, R. G. García, and A. De León Rodríguez, “Estimation of hydrogen production in genetically modified E. coli fermentations using an artificial neural network,” International Journal of Hydrogen Energy, vol. 35, no. 24, pp. 13186–13192, 2010.

J. J. Seppälä, A. Larjo, T. Aho, O. Yli-Harja, M. T. Karp, and V. Santala, “Prospecting hydrogen production of Escherichia coli by metabolic network modeling,” International Journal of Hydrogen Energy, vol. 38, no. 27, pp. 11780–11789, 2013.

Z. Fan, L. Yuan, and R. Chatterjee, “Increased hydrogen production by genetic engineering of Escherichia coli,” PLoS ONE, vol. 4, no. 2, Article ID e4432, 2009.

H. K. Khalil and L. Praly, “High-gain observers in nonlinear feedback control,” International Journal of Robust and Nonlinear Control, vol. 24, no. 6, pp. 993–1015, 2014.