Home | About | Journals | Submit | Contact Us | Français |

**|**PLoS One**|**v.5(9); 2010**|**PMC2940894

Formats

Article sections

Authors

Related links

PLoS One. 2010; 5(9): e12785.

Published online 2010 September 16. doi: 10.1371/journal.pone.0012785

PMCID: PMC2940894

Diego Di Bernardo, Editor^{}

Fondazione Telethon, Italy

Wrote the paper: YJS LB. Conceived the modeling approach and its theoretical foundation: Y-JS. Designed the simulations: Y-JS, LB. Performed the simulations: Y-JS. Analyzed the simulation results: Y-JS, LB.

Received 2010 March 18; Accepted 2010 July 21.

Copyright Shin, Bleris. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

This article has been cited by other articles in PMC.

Systems biology is an interdisciplinary field that aims at understanding complex interactions in cells. Here we demonstrate that linear control theory can provide valuable insight and practical tools for the characterization of complex biological networks. We provide the foundation for such analyses through the study of several case studies including cascade and parallel forms, feedback and feedforward loops. We reproduce experimental results and provide rational analysis of the observed behavior. We demonstrate that methods such as the transfer function (frequency domain) and linear state-space (time domain) can be used to predict reliably the properties and transient behavior of complex network topologies and point to specific design strategies for synthetic networks.

Cells comprise of multiple, heterogeneous subunits that operate in a well-orchestrated manner [1], [2]. Although extremely complex, phenotypes such as cell division and environmental adaptation are the outcome of discrete changes that lead to a deterministic sequence of information transfer and processing within the cell. This information is encoded and transferred via multiple pathways, in different time-scales, and is typically processed in parallel by multi-component networks. Such networks comprise of genes, gene products and small molecules that mutually affect each other or interconvert through biochemical reactions. The number of possible network topologies for a given set of elements is large and it grows exponentially with the number of elements.

Systems biology is an interdisciplinary field that aims at understanding such complex interactions in cells, via the use of a wide spectrum of theoretical and experimental techniques [3]. One of the main thrusts of systems biology is the study of gene networks [4]–[7], via top-down and bottom-up approaches [8]. A top-down approach aims at unraveling the complexity of network dynamics without or with little prior knowledge of the network components and the relationships between them [9]–[12]. On the other hand, a bottom-up approach (closely related to synthetic biology) aims at constructing and studying small-scale biological networks from modular components [13]–[17]. Gene networks are inherently stochastic [18], [19], which renders any modeling effort nontrivial. Furthermore, as the size of a network increases, it becomes increasingly difficult to predict its dynamic behavior. Another characteristic feature is the existence of nonlinearities in biological networks, which further complicates any modeling effort.

Towards analyzing the dynamic behavior of gene networks, a range of mathematical and computational modeling methods have been developed, including Boolean networks, Petri nets, state-charts, ordinary differential equations, and stochastic simulation algorithms [4], [20]–[30]. These approaches can be further organized into two larger categories: logical and continuous models. Logical models deal with the logical sequence of events while continuous models describe the dynamics that depend on molecular concentrations and time.

In this paper, we develop and use linear models [31]–[34], which are shown to capture the dynamics of gene networks in an intuitive and efficient way. We argue that tools of linear control theory, including transfer function (frequency domain) and linear state-space (time domain) methods can be exceptionally practical for systems and synthetic biologists towards unraveling the properties of gene networks and engineering novel systems. We provide several examples of the application of the transfer function method for the analysis of gene networks, starting with network motifs [5], the basic building blocks of gene networks. The transfer function method is sometimes regarded as “classical” in control theory. The state-space or “modern” approach describes a system as a set of input, output, and state variables related by first-order differential or difference equations. One of the advantages of using the state-space method is that it can be used to model multiple-input multiple-output (MIMO) gene networks in a compact manner, utilizing vectors and matrices. Furthermore, we show that the linear state-space approach enables us to utilize a spectrum of tools available for optimal/robust estimation/control for gene network modeling [35]–[37]. As an example, we illustrate that the Kalman filter, one of the well-established optimal estimation tools, can be applied for stochastic modeling of a simple two-gene network. Finally, using a six-node gene network, we demonstrate that our linear approaches can reduce the modeling complexity and provide rapid insight about its dynamic behavior, as compared to conventional non-linear modeling approaches.

The manuscript is organized as follows. We commence our analysis with modeling a simple gene regulation case and subsequently provide the details behind the proposed linearization scheme. Using these results we present the transfer function method for gene network modeling and exemplify the methodology for collapsing cascade and parallel forms to transfer functions. We then provide six case studies that cover a wide range of systems and synthetic biology problems including: cascaded simple regulations, synthetic gene oscillators, effect of the basal production rate, four cascaded simple regulation loops, and finally interconnected feedforward loops. Subsequently, we present the linear state-space method for gene network modeling and we outline its use for the analysis of cascaded simple regulations. Using the linear state-space formulation, we conclude this manuscript with two additional case studies: the use of optimal estimation in gene network measurements and the analysis of the six-node gene network.

A two-gene network (or simple regulation) can be regarded as a fundamental unit that serves as a basic building block for constructing elaborate networks. In simple regulation, one gene (*Y _{gene}*) can be activated by another gene (

(1)

where *y*(*t*) stands for the concentration of *Y _{protein}*. If we assume that

(2)

where *T*_{1/2} is the response time, the time to reach one half of the steady-state *Y _{protein}* concentration. The response times of bacteria are approximately 30 min to a few hours, while those of eukaryotic cells can be longer [13]. Equation 1 can now be restated by substituting

(3)

For the analysis of the function *F(t)*, *w*e need to consider a number of additional factors. As stated earlier, *X _{protein}* must be converted to

(4)

where *x*(*t*) stands for the concentration of total *X _{protein}* that includes both inactive and active forms. It is the maximal level of

(5)

where *x*(*t*) is the maximal level of the *x**(*t*) production that is reached when *s _{x}*(

As a generalization, we can assume that *F*(*t*) is the sum of a basal promoter production rate *F _{0}* and the Hill function:

(6)

*F _{max}* is the maximal level of the

(7)

where *F _{max}* is the maximal level of the

A non-linear model can be linearized around an equilibrium point using a Jacobian matrix, providing insight about the network behavior near that point [38]. For systems that involve Hill functions, there have been various linearization approaches that can be used to study the dynamics outside of the equilibrium [13], [39]–[41]. These methods enable us to convert a non-linear model into a linear model and use linear tools for studying network dynamics over a wide range of parameter and variable values. Two such approaches widely discussed in the literature [13], [39]–[41] are the two-section piecewise linearization (step function or logic approximation) and the three-section piecewise linearization (**Figure 2**). In the figure, *F _{x}*(

Our method is a variation of a two-section piecewise linearization. In the step function approach (yellow in **Figure 2**), *F _{x}*(

The first step towards applying linear control theory tools is to reduce (3) to a linear equation, in other words express *dy*(*t*)/*dt* as a linear combination of x(t) and y(t). The concentration *x**(*t*) of (4) can be expressed as the product of *x*(*t*) and a function *F _{1}*(

(8)

The Hill function relationship between *F _{x}*(

For linearization, we will consider two cases: a case when *x**(*t*) is smaller than *2K _{2A}* and another when it is greater than

(9)

Now, denoting the slope as *F*_{2}, (6) can be expressed as:

(10)

Substituting *x**(*t*) in (10) by *x*(*t*)·*F*_{1}(*t*) from (8), (10) becomes:

(11)

Denoting *F*_{2}·*F*_{1}(*t*) as *f*(t), (11) becomes:

(12)

Substituting *F*(*t*) by *f*(*t*)*x*(*t*), (1) can now be expressed as:

(13)

where dy(t)/dt is shown as a linear combination of x(t) and y(t).

For the *x*(t)* values greater than *2K _{2A}*,

(14)

which again shows that dy(t)/dt is a linear combination of *F _{max}* and y(t). (13) and (14) can be shown together as:

(15)

Similarly, for the case when *s _{x}*(

(16)

When *x*(t)* acts as a repressor, based on the linearization scheme shown in (9–12), the corresponding Hill function (7) can be expressed as:

(17)

From now on, we will consider only the non-trivial cases when *x*(t)* is between zero and 2*K _{2A}*.

A transfer function can be derived from a linear, time-invariant differential equation using Laplace transform [42]. It is a mathematical representation of the relationship between the input and output in the frequency domain. Towards analyzing a gene network, the transfer function method can provide useful information about the behavior and stability of a system. Furthermore, using block diagrams, the transfer function method can represent large complex structures in a simple and intuitive way. Notably, the stochastic nature of gene networks has also been analyzed in the frequency domain [43]–[45].

Starting from equation (8), if we assume that the concentration of the signal *s _{x}*(

(18)

This time-invariance is required for the application of the transfer function method [42]. Using the Laplace transform, (18) can be expressed as:

(19)

where *F(s)*, *X(s)*, and *Y(s)* represent the Laplace transform of *F*(*t*), *x*(*t*), and *y*(*t*). Once the Laplace transform is evaluated, computing in Laplace domain is algebraic and the complexity of solving differential equation is eliminated [42][41][41][42][43][48]. Assuming the initial concentration of *Y _{protein}*,

(20)

Note that the constant *D* plays a critical role in *G*(*s*). It implies that the critical factor that characterizes this system is the dilution/degradation constant. Also in this case the transfer function *G*(*s*) is not affected by the input *F _{0}*. However, as the network becomes more complicated (e.g., with positive and negative feedbacks),

(21)

In case *x(t)* is not changing (has reached its steady state) and set as a constant value *X*, *F·x(t)* becomes *FX* and the resulting step response is:

(22)

where *y(t)* will approach the steady state value as time *t* goes to infinity, and the constant *D* in *e ^{−Dt}* determines how fast the steady state is reached. The step response in Laplace domain can be expressed as:

(23)

Solving a differential equation or evaluating the inverse Laplace transform enables us to evaluate the output response of a system. However, these techniques can be laborious and time-consuming. The use of poles and zeros, fundamental to the analysis and design of control systems, is a technique that can simplify the evaluation. Given the first-order transfer function *G(s)* in (20), a pole, the value for *s* which makes the denominator of the transfer function equal to zero, exists at *s*=−*D*. Since *D* (the degradation/dilution factor) is greater than zero in a real biological system, −*D* can be assumed to be a negative number. This indicates that the system is stable because stable systems do not have any positive poles [42]. For a system to be stable, its natural response must decay to zero, as time approaches infinity, or oscillate. There are various specifications (time constant, rise time, and settling time, etc) that may reveal useful information about the network behavior (the reader is referred to [42]). 1/*D*, called the time constant of the response, is the time for the step response to rise 63% of its final or steady state value. Rise time, defined as the time for the system response to go from 0.1 to 0.9 of its steady state value, is equal to 2.2/*D*. Settling time is the time for the response to reach and stay within 2% of its final value and it is equivalent to 4/*D*.

Based on the block diagram representation of simple regulation shown in **Figure 3**, we can now decompose any complex gene network using four interconnection topologies: cascade forms, parallel forms, feedback loops, and feedforward loops. These topologies are often intermingled with one another as shown in later examples. In this section, we describe the cascade and parallel forms. **Figure 4(a)** shows an example of cascaded simple regulation blocks. The first output is the product of *F(s)* and *G _{Y}(s)*. It is also the input for the second simple regulation. The second or final output is the sum of F(s)G

**Figure 5(a)** shows the block diagram of a feedback loop. Compared to **Figure 3**, there is an additional feedback element *H(s)*. The block diagram can be expressed in the Laplace domain as:

(24)

where *G(s)′* is the simplified equivalent transfer function.

Autoregulation is the simplest form of feedback loops and consists of one simple regulation and one feedback loop (**Figure 5(b)**). In simple regulation, *x**(*t*) is the only transcription factor that determines the production of *y(t)*. On the other hand, in autoregulation, both *x*(t)* and *y*(t)*, the active form of *y(t)*, can affect the *y(t)* production. The *y(t)* production function due to *y*(t)* or *f _{y}(t)* requires two Hill functions as in the case of

(25)

The production function *F _{y}*(

(26)

Applying the linearization scheme used in (9–12), we can rewrite (26) as:

(27)

*h(t)* in (27) is equivalent to *f(t)* in (12). As in (18), if we assume that the signal *s _{y}*(

Negative autoregulation (NAR) occurs when a transcription factor represses the transcription of its own gene (negative feedback). It has been demonstrated that NAR speeds up the response time of the gene expression, decreases the steady state value, and reduces cell-cell variation in protein levels [13], [46]. Using (20) and (24), the transfer function G′(s) can be derived as:

(28)

Note that we have a new degradation/dilution constant *D′*. Since both *D* and *H* have positive values in biological systems, *D′* is also positive, indicating that the system is stable based on the pole-zero analysis described previously (**Figure 5(c)**). Furthermore, the time constant of the response (*1/D′*) is greater than *1/D* since *H* is positive. This explains the experimental observation of the decrease of the response time using NAR (**Figure 5(c)** and **(d)**).

Positive autoregulation (PAR) occurs when a transcription factor controls its own protein production rate. In contrast to NAR, the response time is extended and the steady state value and cell-cell variation are increased [13]. Using (20) and (24), the transfer function of PAR is:

(29)

Equation 29 illustrates that, compared to simple regulation, the time constant (*1/D″*) is increased because *D*″ is equal to *D* subtracted by *H*, a positive number. This explains the increase of the response time (**Figure 5(c)** and **(d)**). Additionally, for the system to be stable, – *D″* must be negative (or *H* must be smaller than *D*). In other words, if the positive feedback is “too strong” (*H* is greater than *D*) then the system may become unstable.

**Figure 6** shows the block diagram of a negative feedback loop that consists of two cascaded simple regulations. *X _{gene}* activates

(30)

Substituting *G _{Y}(s)* and

(31)

Now we have a second-order transfer function *G(s)*. Whereas varying the first-order transfer function parameter *D* changes only the response time (1/*D*), changing the parameters (*F _{YZ}*,

(32)

The natural frequency *ω _{n}* of a second-order system is defined as the frequency of an undamped oscillation [42].

(33)

Without damping, the poles would be on the imaginary axis as shown in **Figure 7(b)**. Equation (33) shows that the natural frequency can be tuned by varying the strength of both *F _{YZ} and F_{ZY}*.

For a damped system, the damping ratio *ζ* is defined as [42]:

(34)

Various step responses of a second-order system with respect to the damping ratio *ζ* are shown in **Figure 7**. The damped natural frequency *ω _{n}* is defined as [42]:

(35)

The poles of an undamped oscillating system lie on the imaginary axis, and thus the system is said to be marginally stable [42]. Also, note that all the poles in **Figure 7** lie on the left half of the plane except for the ones that belong to an undamped oscillation. For an underdamped oscillation, additional specifications can be defined and used to predict the response in detail, including the peak time, percent overshoot, settling time, and rise time. Using a numerical approach, rise time can be computed even though there is no precise analytical solution [42]. All these specifications can be expressed in terms of *F _{YZ}*,

Building tunable synthetic gene oscillators has been a central area of focus for systems and synthetic biologists [47]–[49]. An oscillating system is basically an undamped second-order transfer function (**Figure 7(b)**). For a system to have no damping, *a* (being equal to *D _{Y}* +

A tunable synthetic gene oscillator is shown in **Figure 8** [47]. The genes araC (denoted as 1) and lacI (denoted as 2) have identical hybrid promoters that can be activated by AraC in the presence of arabinose and repressed by LacI in the absence of IPTG. As shown in the figure, araC has a positive autoregulation while LacI a negative autoregulation. *F _{1}* and

(36)

According to (36), *a* (=*D _{1}* +

(37)

Equation 37 shows that both the parameters influenced by arabinose (*F _{1}* and

Using the transfer function method, we can identify a way of changing such a non-monotonic relationship into a monotonic one. For example, if the negative autoregulation for *lacl* is removed from the network (by removing the repression operator site from the promoter), *D _{1}F_{3}* and −

(38)

Now, *F _{4}* belongs to a positive term, meaning that as the IPTG concentration is raised the oscillating frequency also increases, and vice versa, indicating that IPTG can now cause a monotonic behavior.

In order for an undamped system (*a*=0) to be marginally stable, *b* in (37) must be greater than zero so that the two poles lie separately on the imaginary axis [42]:

(39)

There are two negative terms (– *D _{2}F_{1}* and –

In this section, we will demonstrate how the effect of the basal production rate *f _{0}*(

(40)

Note that *E*(*s*) is the error (difference) between F(s) and *F _{ZY}Z*(

(41)

where the first term can be regarded as a transfer function relating *E*(*s*) to *F*(*s*) and the second term relating *E*(*s*) to *F _{0,Z}*(

Applying the final value theorem to Eq. (41), we obtain:

(42)

As we assume the basal rate is constant *F _{0,Z}*,

(43)

Equation (43) can provide important insight about reducing the steady-state error due to basal activity by tuning the strength of *F _{YZ}* or

A second order system with a negative feedback loop described earlier could be efficiently analyzed using a pole/zero plot and various formulas. Even though we do not have such formulas for higher order systems, pole/zero plots can still be useful for predicting their behavior.

**Figure 10** shows the block diagram of a fourth order system with various feedbacks. It is a synthetic network called IRMA (*In vivo* Reverse-engineering and Modeling techniques Assessment) that consists of four genes (*CBF1*, *GAL4*, *SWI5*, and *ASH1*) in *Saccharomyces Cerevisiae* [50]. It has been demonstrated both computationally (using a non-linear model) and experimentally that the network could be turned into an oscillator by changing various parameters (Michaelis-Menten coefficient, Hill coefficient, etc) that control the strength of the interactions among genes.

An equivalent model was built using the transfer function method, as shown in **Figure 10(a)**, where the strength of the interactions was determined by changing the values of the constants (F_{12}, F_{23}, F_{31}, F_{33}, F_{34}, and F_{41}) shown in the figure. The numbers 1, 2, 3, and 4 in subscript represent the four genes mentioned earlier, respectively. As described in the paper [50], we were able to demonstrate that increasing the values of F_{12}, F_{33}, F_{34}, and F_{41} contributes to generating an oscillatory behavior. **Figure 10(b)** shows that as we increase the values of those constants, the dominant poles (two poles that have the least real values) cross the imaginary axis. When the dominant poles are exactly on the axis, the system exhibits an undamped oscillation, similar to the behavior of a second order system shown in **Figure 7** (see Movie S1 and File S1).

Feedforward loops (FFL) are network motifs [13] that combine into more complex and larger networks (e.g. in *Bacillus Subtilis* [51]). The coherent type-1 feedforward loop (C1-FFL) and incoherent type-1 feedforward loop (Ic1-FFL) are the most abundant FFL types [52]. In this section, we will describe how transfer function method can be used to model a network that consists of interconnected feedforward loops. **Figure 11(a)** shows a simplified schematic diagram of the network. Two sets of C1-FFL and Ic1-FFL in parallel are connected in cascade. Ic1-FFLs generate pulses of C1 and C2, and C1-FFLs create delays in C2 and C3 expression. **Figure 11(b)** shows a schematic diagram that illustrates the sequential expression of C1, C2, and C3. **Figure 11(c)** shows the block diagram of the network. Using the transfer function method, the sequential expression of the three genes can be simulated as shown in **Figure 11(d)**. Note that if the values of the protein concentrations are normalized, the expression pattern will be equivalent to the one shown in **Figure 11(b)**.

A state is a complete summary of the status of a system at a particular point in time, and it is described by the values of a set of state variables [53]. Based on the linearization scheme presented previously, a linear state-space method can be applied to model gene networks. There are many benefits in using the linear state-space approach. First, it allows time-varying systems, meaning that the models can include time-varying signals. This was not possible in the case of transfer function-based models described earlier. Secondly, large complex networks that are multiple-input multiple-output (MIMO) systems can be represented in a compact way, using vectors and matrices. Thirdly, the dynamic behavior of a system can be understood using the eigenvalues and eigenvectors. The eigenvalues are equivalent to the poles discussed in the transfer function method. Finally, we can build a stochastic linear state-space model that enables us to utilize a spectrum of tools available for optimal/robust estimation/control for gene network modeling [35]–[37]. As an example, we illustrate that the Kalman filter, one of the well-established optimal estimation algorithms in science and engineering, can be applied towards modeling a simple two-gene network.

In the linear model of simple regulation (linear ODE form shown in (13)), the system involves two proteins, *x*(*t*) and *y*(*t*), with units of concentration per cell. If we consider gene *y* as our system of interest, then and the state vector **s** can be expressed as:

(44)

The state-space model can be represented as [53]:

(45)

where ** u** is the input vector. Note that we do not have to assume that

The state-space method is analogous to the transfer function method in many ways. Therefore, in this section we present only the case of a negative feedback loop involving two cascaded simple regulations. The network shown in **Figure 6** can be written as a state-space model:

(46)

where *f _{ab}*(

(47)

Since the eigenvalues are equivalent to the poles, we can immediately predict that *y* and *z* will exhibit underdamped oscillations and the system is stable. The solution can be written as:

(48)

where *C*_{1} and *C*_{2} are determined by the initial values of *y* and *z*. Depending on the initial values, the time-dependent values of *y* and *z* can follow an infinite number of paths as shown in the vector field (**Figure 12**). However, regardless of the initial values, the figure shows that all the trajectories will eventually converge into an equilibrium point. This equilibrium point can also be computed by solving *dy*(*t*)/*dt*=0 and *dz*(*t*)/*dt*=0 in (46). Note that the damped oscillatory behaviors that were predicted using the eigenvalues are observed in both *Y* and *Z* expressions (**Figure 12**).

Stochastic optimal estimation [18], [19], [43], [45] and control methods aim to determine the best strategy for estimating or controlling a dynamic system in the presence of uncertainty [35]. This objective is common to many fields, including engineering, science, and economics [35]. Assume that, in the case of simple regulation x→y, we are given noisy time-series fluorescence measurement data of *y* protein [55]–[57]. Would it be possible to “optimally” estimate protein *x* which is also fluctuating due to noise? In this section, we will illustrate that the Kalman filter, a well-established optimal estimation algorithm, can be applied for stochastic modeling of a simple regulation (two-gene network). Before we discuss the stochastic model of simple regulation, it is important to note that there are two different noise contributions in gene networks [18], [58], [59]. “Intrinsic noise” is generated by the inherent stochasticity of biochemical processes, such as transcription and translation, which are directly related to the expression of a specific gene. On the other hand, the environment including other cellular components (mitochondria, microtubules, etc) that indirectly influence the expression of *y* gene, contribute to the “extrinsic noise”.

First, the least squares estimation is used to estimate the constant concentration of protein *x* (without noise), given noisy measurement data of protein *y*. Using the Euler's method, equation (13) can be expressed as:

(49)

where *h* is the step size. For simplicity, we will assume *f*(*t*) and *d*(*t*) are constant (*F* and *D*) and the basal protein production rate and initial value of *y* are zero (*F*_{0,n}=*y*_{1}=0):

(50)

which can be compactly shown as ** y**=

(51)

where *e** ^{y}* is assumed to be a Gaussian zero-mean white noise vector.

(52)

Note that in this formulation protein x is noise-free and thus has no stochastic influence on y. Using the least squares method, the estimated ** x** can be calculated as:

(53)

Now, we will use the Kalman filter for estimating noisy ** x** from noisy

(54)

**Figure 13(a)** shows ** x**, which is the sum of a constant value (200 molecules/cell) vector and Gaussian zero-mean white noise vector with a standard deviation of 10. The noise is assumed to have both intrinsic and extrinsic components.

Here, we will illustrate the utility of the state-space method for the analysis of a 6-node network (**Figure 14(a)**). There are a few dynamic features that we may intuitively extract from the figure based on our previous discussions. For example, there is a coherent type-1 feedforward loop (C1-FFL) embedded in the network (that consists of A, B, and C), which suggests that gene C is expressed after the expression of gene B (there is a delay in the C expression).

The 6-node gene network can be modeled using the following state-space formulation:

(55)

where *f _{AB}*(

We demonstrated that well-established tools of linear control theory can be used to model gene networks and explain experimental observations in a highly intuitive way. Methods such as the transfer function (frequency domain) and linear state-space (time domain) were applied to reveal inherent characteristics and predict the dynamic behavior of various gene network topologies, including cascade/parallel forms, feedback loops, and feedforward loops. Additionally, we showed that well-established optimal estimation tools, such as the Kalman filter, can be used in the context of gene network modeling in the presence of noise.

While we assumed that multiple transcription factors act on a single gene in an additive/subtractive way, in biological systems positive or negative cooperativity, and mutual exclusion can also be observed. Combinational logic is one approach that can be used for determining the net effect of multiple inputs at such junctions [13], [30]. We expect that hybrid models that combine both continuous and logical approaches may provide additional insight.

For simulations, the value of degradation/dilution constant (D) was 0.01/min and *f*(*t*) was set as a constant (*F*) that ranged from 0.01 to 0.1/min. The input *x*(*t*) had a fixed value between 100 and 1000 molecules/cell and the typical basal protein production rate was set as 0.1 molecules/min.

For the estimation of ** x**, a new piece of information y

(56)

From (56), it can be stated that:

(57)

A linear recursive estimator can be written as [36]:

(58)

where *K*_{N} is a gain matrix.

Suppose we have the following linear discrete-time system [36]:

(63)

where is a known input and is Gaussian zero-mean white noise. The covariance matrix of *e*^{x} can be shown as:

(64)

If we take the expected value of both sides of (63):

(65)

which shows the propagation of the mean of *x** _{N}* [36]. The propagation of the covariance can be shown as [36]:

(66)

If we have all of the measurements up to and including time N+1 available for our estimate of , then we can form *a posteriori* estimate denoted as . If we have all of the measurements before (not including) time N+1 available for our estimate of , then we can form *a priori* estimate denoted as _{.} Similarly, denotes the covariance of and denotes the covariance of . Assuming that the initial state is given:

(67)

Using (65) that describes the propagation of the mean of , we obtain:

(68)

The reasoning can be extended as following:

(69)

Assuming , the covariance of the initial estimate , is given and using (66):

(70)

The reasoning can be extended as following:

(71)

Now we have the time-update equations for and *P* and we need the measurement-update equations for and *P*. Utilizing the recursive least squares development (56–62), the discrete-time Kalman filter for simple regulation can be summarized as:

(72)

LabVIEW VI file for Case Study V.

(0.06 MB TXT)

Click here for additional data file.^{(61K, txt)}

**Competing Interests: **The authors have declared that no competing interests exist.

**Funding: **These authors have no support or funding to report.

1. Nelson DL, Cox MM. Lehninger Principles of Biochemistry. New York: W H Freeman and Company; 2005.

2. Alberts B, Johnson A, Lewis J, Raff M, Roberts K, et al. Molecular Biology of the Cell. New York: Garland Science; 2002.

3. Kitano H. Systems biology: a brief overview. Science. 2002;295:1662–1664. [PubMed]

4. Hasty J, McMillen D, Isaacs F, Collins JJ. Computational studies of gene regulatory networks: in numero molecular biology. Nat Rev Genet. 2001;2:268–279. [PubMed]

5. Milo R, Shen-Orr S, Itzkovitz S, Kashtan N, Chklovskii D, et al. Network motifs: simple building blocks of complex networks. Science. 2002;298:824–827. [PubMed]

6. Iglesias PA, Ingalls BP. Control Theory and Systems Biology. Cambridge: The MIT Press; 2009.

7. Alon U. Network motifs: theory and experimental approaches. Nat Rev Genet. 2007;8:450–461. [PubMed]

8. Cuccato G, Della Gatta G, di Bernardo D. Systems and Synthetic biology: tackling genetic networks and complex diseases. Heredity. 2009;102:527–532. [PubMed]

9. Schena M. Quantitative monitoring of gene expression patterns with a complementary DNA microarray. Science. 1995;270:467. [PubMed]

10. Kohane IS, Kho AK, Butte AJ. Microarrays for an Integrative Genomics. Cambridge: The MIT Press; 2005.

11. Sprinzak D, Elowitz MB. Reconstruction of genetic circuits. Nature. 2005;438:443–448. [PubMed]

12. DasGupta B, Vera-Licona P, Sontag E. Reverse Engineering of Molecular Networks from a Common Combinatorial Approach. In: Elloumi M, Zomaya AY, editors. Algorithms in Computational Molecular Biology: Techniques, Approaches and Applications. Hoboken: John Wiley & Sons; 2010.

13. Alon U. An Introduction to Systems Biology: Design Principles of Biological Circuits. London: CRC; 2006.

14. Guido NJ, Wang X, Adalsteinsson D, McMillen D, Hasty J, et al. A bottom-up approach to gene regulation. Nature. 2006;439:856–860. [PubMed]

15. Andrianantoandro E, Basu S, Karig DK, Weiss R. Synthetic biology: new engineering rules for an emerging discipline. Mol Syst Biol. 2006;2:2006.0028. [PMC free article] [PubMed]

16. Bleris LG, Xie Z, Adadey A, Sontag E, Benenson Y. Transcriptional feedforward motif as stable expression unit. 2009. Systems Biology of Human Disease Conference Proceeding. Boston.

17. Lu TK, Khalil AS, Collins JJ. Next-generation synthetic gene networks. Nat Biotechnol. 2009;27:1139–1150. [PMC free article] [PubMed]

18. Dunlop MJ, Cox RS, 3rd, Levine JH, Murray RM, Elowitz MB. Regulatory activity revealed by dynamic correlations in gene expression noise. Nat Genet. 2008;40:1493–1498. [PMC free article] [PubMed]

19. Raj A, van Oudenaarden A. Nature, nurture, or chance: stochastic gene expression and its consequences. Cell. 2008;135:216–226. [PubMed]

20. Thomas R. Boolean formalization of genetic control circuits. J Theor Biol. 1973;42:563–585. [PubMed]

21. Smolen P, Baxter DA, Byrne JH. Modeling transcriptional control in gene networks–methods, recent results, and future directions. Bull Math Biol. 2000;62:247–292. [PubMed]

22. Bower JM, Bolouri H. Computational Modeling of Genetic and Biochemical Networks. Cambridge: MIT Press; 2001.

23. Shmulevich I, Dougherty ER, Zhang W. Gene perturbation and intervention in probabilistic Boolean networks. Bioinformatics. 2002;18:1319–1331. [PubMed]

24. Goutsias J, Lee NH. Computational and experimental approaches for modeling gene regulatory networks. Curr Pharm Des. 2007;13:1415–1436. [PubMed]

25. Steggles LJ, Banks R, Shaw O, Wipat A. Qualitatively modelling and analysing genetic regulatory networks: a Petri net approach. Bioinformatics. 2007;23:336–343. [PubMed]

26. Schilstra Bio-logic: Gene expression and the laws of combinatorial logic. Artif Life. 2008;14:121. [PubMed]

27. Karlebach G, Shamir R. Modelling and analysis of gene regulatory networks. Nat Rev Mol Cell Biol. 2008;9:770–780. [PubMed]

28. Garg A, Di Cara A, Xenarios I, Mendoza L, De Micheli G. Synchronous versus asynchronous modeling of gene regulatory networks. Bioinformatics. 2008;24:1917–1925. [PubMed]

29. Bolouri H. Computational Modeling of Gene Regulatory Networks. London: Imperial College Press; 2008.

30. Shin YJ, Nourani M. Statecharts for gene network modeling. PLoS One. 2010;5:e9376. [PMC free article] [PubMed]

31. Anderson BDO, Moore JB. Optimal Control: Linear Quadratic Methods. Mineola: Dover; 1990.

32. Lathi BP. Linear Systems and Signals. Oxford: Oxford University Press; 2002.

33. Schweickhardt T, Allgöwer F. Quantitative nonlinearity assessment – An introduction to nonlinearity measurement. Amsterdam: Elsevier Science; 2004.

34. Strang G. Computational Science and Engineering. Wellesley: Wellesley-Cambridge Press; 2007.

35. Stengel RF. Optimal Control and Estimation. Mineola: Dover; 1994.

36. Simon D. Optimal State Estimation. New York: John Wiley & Sons; 2006.

37. Skogestad S, Postlethwaite I. Multivariable Feedback Control. New York: John Wiley & Sons; 2005.

38. Iglesias PA, Ingalls BP. Control Theory and Systems Biology. Cambridge: MIT Press; 2010.

39. De Jong H, Gouze JL, Hernandez C, Page M, Sari T, et al. Qualitative simulation of genetic regulatory networks using piecewise-linear models. Bull Math Biol. 2004;66:301–340. [PubMed]

40. Casey R, de Jong H, Gouze JL. Piecewise-linear models of genetic regulatory networks: equilibria and their stability. J Math Biol. 2006;52:27–56. [PubMed]

41. Polynikis A, Hogan SJ, di Bernardo M. Comparing different ODE modelling approaches for gene regulatory networks. J Theor Biol. 2009;261:511–530. [PubMed]

42. Nise NS. Control Systems Engineering. New York: John Wiley & Sons; 2000.

43. Simpson ML, Cox CD, Sayler GS. Frequency domain analysis of noise in autoregulated gene circuits. Proc Natl Acad Sci U S A. 2003;100:4551–4556. [PubMed]

44. Cox CD, McCollum JM, Austin DW, Allen MS, Dar RD, et al. Frequency domain analysis of noise in simple gene circuits. Chaos. 2006;16:026102. [PubMed]

45. Austin DW, Allen MS, McCollum JM, Dar RD, Wilgus JR, et al. Gene network shaping of inherent noise spectra. Nature. 2006;439:608–611. [PubMed]

46. Becskei A, Serrano L. Engineering stability in gene networks by autoregulation. Nature. 2000;405:590–593. [PubMed]

47. Stricker J, Cookson S, Bennett MR, Mather WH, Tsimring LS, et al. A fast, robust and tunable synthetic gene oscillator. Nature. 2008;456:516–519. [PubMed]

48. Tsai TY, Choi YS, Ma W, Pomerening JR, Tang C, et al. Robust, tunable biological oscillations from interlinked positive and negative feedback loops. Science. 2008;321:126–129. [PMC free article] [PubMed]

49. Tigges M, Marquez-Lago TT, Stelling J, Fussenegger M. A tunable synthetic mammalian oscillator. Nature. 2009;457:309–312. [PubMed]

50. Marucci L, Barton DA, Cantone I, Ricci MA, Cosma MP, et al. How to turn a genetic circuit into a synthetic tunable oscillator, or a bistable switch. PLoS One. 2009;4:e8083. [PMC free article] [PubMed]

51. Eichenberger P, Fujita M, Jensen ST, Conlon EM, Rudner DZ, et al. The program of gene transcription for a single differentiating cell type during sporulation in Bacillus subtilis. PLoS Biol. 2004;2:e328. [PMC free article] [PubMed]

52. Mangan S, Alon U. Structure and function of the feed-forward loop network motif. Proc Natl Acad Sci U S A. 2003;100:11980–11985. [PubMed]

53. Brogan WL. Modern Control Theory. Englewood Cliffs, New Jersey: Prentice-Hall; 1991.

54. Luenberger DG. Introduction to Dynamic Systems. New York: John Wiley & Sons; 1979.

55. Suel GM, Garcia-Ojalvo J, Liberman LM, Elowitz MB. An excitable gene regulatory circuit induces transient cellular differentiation. Nature. 2006;440:545–550. [PubMed]

56. Dunlop MJ, Cox RS, 3rd, Levine JH, Murray RM, Elowitz MB. Regulatory activity revealed by dynamic correlations in gene expression noise. Nat Genet. 2008;40:1493–1498. [PMC free article] [PubMed]

57. Locke JC, Elowitz MB. Using movies to analyse gene circuit dynamics in single cells. Nat Rev Microbiol. 2009;7:383–392. [PMC free article] [PubMed]

58. Raser JM, O'Shea EK. Noise in gene expression: origins, consequences, and control. Science. 2005;309:2010–2013. [PMC free article] [PubMed]

59. Swain PS, Elowitz MB, Siggia ED. Intrinsic and extrinsic contributions to stochasticity in gene expression. Proc Natl Acad Sci U S A. 2002;99:12795–12800. [PubMed]

Articles from PLoS ONE are provided here courtesy of **Public Library of Science**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |