Yang et al. 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.
External Control of the GAL Network in S. cerevisiae: A View from Control Theory
1Department of Mechanical, Aerospace and Biomedical Engineering, The University of Tennessee, Knoxville, Tennessee, United States of America
2Vanderbilt Institute for Integrative Biosystems Research and Education, Departments of Biomedical Engineering, Molecular Physiology & Biophysics, and Physics & Astronomy, Vanderbilt University, Nashville, Tennessee, United States of America
Diego Di Bernardo, Editor
Fondazione Telethon, Italy
Received September 15, 2010; Accepted March 31, 2011.
While there is a vast literature on the control systems that cells utilize to regulate their own state, there is little published work on the formal application of control theory to the external regulation of cellular functions. This paper chooses the GAL network in S. cerevisiae as a well understood benchmark example to demonstrate how control theory can be employed to regulate intracellular mRNA levels via extracellular galactose. Based on a mathematical model reduced from the GAL network, we have demonstrated that a galactose dose necessary to drive and maintain the desired GAL genes' mRNA levels can be calculated in an analytic form. And thus, a proportional feedback control can be designed to precisely regulate the level of mRNA. The benefits of the proposed feedback control are extensively investigated in terms of stability and parameter sensitivity. This paper demonstrates that feedback control can both significantly accelerate the process to precisely regulate mRNA levels and enhance the robustness of the overall cellular control system.
The complexity of nested feedback control loops associated with complex biological systems creates many challenges to understanding biological systems. At the same time, the need for precise control of gene expression is increasing in fields such as genetics, cell and molecular engineering, and in disease treatment via gene expression control. A straightforward application for control of gene expression might be to generate a maximum output in bioreactor-based systems. The outputs from this type of application would typically by specific proteins of interest for biopharmaceuticals. The approach can be further applied to create a system in which the control of a gene network may be coupled with the expression of a recombinant protein.
On the other hand, systems biologists have built numerous mathematical models for gene and protein networks that cells utilize to control themselves 
. The applications of such models have gone beyond simple qualitative understanding of network dynamics. Some models have been used in synthetic biology to affect metabolism and eventually control the biological processes toward desired outputs 
. The full advantages of control theory, however, have not been realized, and there is little work on external control of cellular functions such as gene expression.
From a control theory perspective, control techniques used for complex engineering systems should in principle be applicable to the regulation of cellular systems. There are in general two potential challenges. First, it is not easy to find simple mathematical models for cellular system control. Most existing cellular mathematical models are not formulated in a standard control system framework. Typically, they lack global feedback loops from outputs to dynamically adjust the control dosage. Second, model parameters are usually very sensitive due to large uncertainties. The resulting control algorithm often becomes fragile, and the control system can remain stable only under small disturbances. Open-loop constant controls – for example, the gene knockdown or overexpression technique for yeast performed without adjusting dosages in real-time in response to the observed changes – are vulnerable to disturbances and individual heterogeneity.
It is desirable to have model-based, fine-tuned external control approaches to precisely regulate a reverse-engineered target. Unfortunately, achieving a reliable feedback control remains a challenge for cellular systems. Recently, Cantone et al.
developed the first reverse-engineering benchmark system 
, which is based on the GAL network using real-time image feedback. Control engineers have begun to apply closed-loop control ideas based on this benchmark 
In this paper we have employed model reduction techniques to simplify the complex GAL network, and further demonstrate that the simplified model is favorable for control system design. We also present the benefit of closed-loop system design. Our goal is to demonstrate the usefulness of control theory for precise external regulation of cellular systems and to open discussions in the field about possible future benefits of more advanced theoretical control.
network in yeast, Saccharomyces cerevisiae
, is one of the most studied gene networks, and hence is well suited to serve as a demonstration, benchmark system. Galactose utilization in the network can be broken up into four distinct stages based on proteins present in the cytoplasm. In the first stage GAL
2p, the galactose transporter protein, allows galactose to enter into the cytoplasm 
. Once internalized, GAL
3p binds to the cytoplasmic galactose and is activated. This causes the galactose to bind and inhibit the action of GAL
. With GAL
80p bound, GAL
4p is able to activate the transcription of GAL
3, and GAL
80 genes. In the presence of glucose, GAL
4p is also inhibited, which leads to inhibition of GAL
1 expression 
. Without accounting for the glucose network, the GAL
system is composed of two positive feedback loops, GAL
3p and GAL
2p, and two negative feedback loops, GAL
80p and GAL
1p. The glucose network proteins bind upstream GAL
3, and GAL
4, and suppress transcription through the action of Mig1. The glucose network also prevents the binding of GAL
2p and external galactose and shuts down the galactose transportation. A mathematical model of the galactose network has previously been established by de Atauri et al. 
and Ramsey et al. 
. Bennett et al.
combined a simplified glucose network with the galactose network model 
, and identified the model parameters using experimental data. A diagram of the simplified glucose and galactose networks used in this study is shown in .
Figure 1 The gene regulatory network for galactose utilization (redrawn from Bennett et al. (2008) ).
The key factor in designing a control system for a complex network that can be controlled externally is to design the system with measurable input(s) and output(s). For our study of the GAL
network, the measurable inputs into the system are the concentrations of galactose and glucose. The most common mechanism currently available to monitor the output is a reporter to a gene of interest, normally green fluorescent protein (GFP), which allows researchers to measure the expression level of the gene. In the theoretical control model designed in this study, a GFP reporter coupled to GAL
1 can be used to measure the output from the system. In practice, a spectrofluorometer or a microscope with appropriate excitation and emission filters, possibly coupled to a microfluidic device, could be used for measuring the GFP concentration 
. It would then be possible to develop a control algorithm for the expression of GAL
1 by tuning the concentration of galactose delivered externally to the cells. This model could be expanded to the control of other genes of interest (GAL2
3, and GAL
80) coupled to separate reporters.
Another challenge for control system design is the network scale and model complexity 
. A gene network model is often constructed based on chemical kinetics with tens or hundreds of states and many parameters 
. The most successful model for the yeast cell cycle alone has 61 coupled ordinary differential equations and 141 parameters 
. In such models, many important system characteristics, such as steady states, cannot be solved analytically and have to be estimated using numerical simulation. Advanced control mechanisms directly applied to these systems may waste considerable time on calculation of complicated, yet inessential nonlinear reaction terms. This complexity may lead to fragility of the controlled system, causing the system to collapse due to noise and parameter perturbation. The reduction of model complexity thus becomes extremely important for complex control system design and implementation. One effective tool for model reduction is parameter sensitivity analysis, which elucidates the dependence of system dynamics on the parameters. A small sensitivity measure for a parameter implies that the value of this parameter can be substituted for a wide range of values without altering the system dynamics. Some reaction terms can thus be deemed negligible and result in a reduced model. The control design based on the reduced model can be expected to have a similar performance to that of the original model. This paper uses the GAL
1 mRNA level as a measureable output to be controlled (with the assumption that mRNA level is highly correlated with GAL
1 translation) and investigates the sensitivity between the output as well as model parameters. The resulting reduced model based on sensitivity analysis successfully separates the galactose utilization network and the glucose network, and reduces the GAL
4p subsystem and associated complexes. Furthermore, global sensitivity was also conducted to show how the system changes due to the simultaneous variation of all the parameters over a wide range of values. We have concluded that the global sensitivity analysis is consistent with the results of local sensitivity. In addition, the nonlinearity of the system is largely reduced, while maintaining a small deviation in the system dynamics. After system analysis, the original complex system was simplified to allow for the control system design. The goal of the control is to maintain one of the GAL
mRNAs at a desired level.
In this paper, both a constant open-loop control and a proportional-output feedback closed-loop control are designed based on the reduced model. First, an analytic formula for the precise dosage of the external galactose is given, instead of manually tuning galactose by trial and error. Feedback control is introduced to improve the robustness of the controlled system and enhance the convergence rates. The simulation results show that both controls can achieve the control objective, i.e., maintaining GAL1 mRNA at a desired level. Similar analytic galactose dosages can be achieved for all the other measurable outputs, such as GAL2, GAL3, and GAL80. The difference between the open-loop control of the system with constant control input, and feedback control with time-dependent parameters, is that the feedback control significantly shortens the time required to achieve the steady state. The feedback loop indirectly increases the degradation rate of the internal galactose so that the system more rapidly reaches balance. The feedback control also significantly reduces parameter sensitivity. Feedback control takes advantage of information about the system state to regulate output, and in doing so significantly decreases the local sensitivity of measurements to the system parameters. Hence feedback control can resist much larger parameter perturbations/uncertainties as compared to the constant open-loop control; that is, an open-loop control with constant control input.
Mathematical Model for Yeast GAL Network
network is a good starting point to demonstrate feedback and feed forward external controls for cellular systems. Though simple, it is well understood and is easily manageable to interpret experimental results. On the other hand, it is complex enough to test sophisticated control algorithms. As illustrated in , a mathematical model for the GAL
network has been proposed based on the interactions between proteins and internalized galactose. The figure shows a gene regulatory network for galactose utilization (redrawn from Bennett et al.
). The extracellular inputs are [gal
. The natural output of the network is a group of proteins. The GAL
4 protein g4
binds to upstream activation sites and activates the regulatory genes in the galactose network. The GAL
80 gene inhibits the inducing effects of GAL
4 and thereby provides negative feedback in the system. GAL
3 enhances expression of GAL
4 by binding with internal galactose (Gal), forming a GAL
3-galactose complex, g3c
, that inactivates g80
by binding to it and resulting in a complex g80c
. In addition, the transporter GAL
2 increases the amount of internal galactose, which stimulates the galactose network. We use the GAL
1 mRNA level, m1
, as the measureable output to control the system. According to sensitivity analysis, the interaction loops involving GAL4p dimerization (the dotted lines) can be eliminated from the original model to create a reduced model. In the figure, gi
is the number of the galactose network protein monomers (i
1, 2, 3, 4, 80), mi
is amount of mRNA (i
1, 2, 3, 80), gid
is the number of protein dimers (i
is the number of GAL
3p proteins bound to galactose,
is the number of GAL
4p dimers bound to Gal80p dimers, and
is the number of GAL
80p proteins bound to the Gal3p-galactose complex.
A complete mathematical model of the above gene network, which we term the Original model, can be described with 22 mathematical equations 
. The first five represent the mass-action kinetics of galactose protein monomers, including dimerization, and the interaction with the internal galactose:
Equations (6)–(9) describe mRNA kinetics accounting for the transcription and translation of the galactose genes as well as the degradation of mRNA
Equations (10) – (14) represent the kinetics of protein dimers and the associated complexes:
The metabolic reactions and transport of galactose are represented by a single equation:
Equations (16) – (18) describe a simplified glucose network, including the glucose-mediated enzymatic decay of GAL1
, and [glu
] are the amount of glucose network mRNA, associated protein, and cellular internal glucose, respectively. They compose a simplified glucose network. The inhibitory effect due to products of the glucose network that act on various processes of the galactose network is represented as follows,
are the transport rates of external galactose and external glucose into the cell, respectively. The cooperative fractional saturation function describing the number of upstream activation sites occupied on a promoter, assuming that N
sites exist 
, is given by
The above mathematical model is mainly based on the biochemical reactions occurring throughout the network. This approach of modeling the intermediates in a pathway has been widely used in both chemistry and biology. The variables and parameters used for the model are defined and summarized in and .
Our goal for galactose network control is to regulate the GAL
family mRNA level by manipulating galactose and glucose concentrations. In this paper, we select GAL
1 as the gene of interest. The mRNA level of this target gene can be measured experimentally using a GFP conjugated to GAL
1, which can be measured for intact cells, or by means of a microarray or other assay that requires lysing a small fraction of the cells under observation. In practice, other reporter genes may be used to monitor GAL
1 expression, for example, the red fluorescent protein from the gene dsRed 
. In addition, it is notable that glucose suppresses all GAL
. External glucose must be kept at a very low level, i.e
, or the GAL
network will be inhibited. Once no external glucose is supplied, the glucose network will soon reach its steady state, which means one can focus on the galactose network.
From a control engineering perspective, one key step is to identify inputs (control) and outputs (preferably measurable), so that the system can be put into a control framework for discussion.
, external galactose and glucose levels, are the control variables in this study, while the mRNA levels of GAL
3, and GAL
, are measurable variables. The goal of this research was to design a time-course of
such that one of the mRNAs m1
is maintained at a desired level.
While there are several methods that might be suitable for culturing yeast under external control, microfluidic devices would be ideal for experimental implementation of the cellular control system 
. Fluorescence microscopy could be used to quantify, in real-time, the levels of GFP-labeled Gal
1p from a small, synchronized population of S. cerevisiae,
with the output of a camera or photomultiplier tube serving as the sensor measurement to the control system. Ion mobility-mass spectrometry could also be used to monitor the secreted metabolites and signaling molecules in real-time and intracellular species after cell lysis 
. Microfluidic valves could control the concentrations of various chemicals and serve as the control system outputs.
In the next section, we will conduct a sensitivity analysis to simplify this subsystem. As shown in , the interaction loops involving GAL4p dimerization can be removed. The reduced model lowers the complexity to design the associated controller, allowing analytical implementation of a feasible design. We will show later that the reduced model may cause bounded deviation from the original model in transient response; however, it results in insignificant deviation from the steady states.
Given the steady states of the glucose utilization network, the original model still has 35 parameters. The large number of parameters and significant nonlinearities make it difficult to design an effective controller. As with any biological system, the number of parameters involved in the induction of a specific pathway could be large, with physical, chemical, and biological parameters all affecting the system. In terms of a modeling-based representation of biological systems, the goal of our effort is to determine the parameters that most affect the induction of the system.
A simple example of this concept can be observed in glycolysis. Through the complex intermediates and interactions within this pathway, three regulatory steps are often considered, biologically, the most important rate-limiting steps in the pathway. These irreversible steps, phosphorylation of glucose, phosphorylation of fructose-6-phosphate, and transfer of phosphate to phosphoenolpyruvate and then to ADP, would be the most important control points. The proteins involved in these steps, as well as the mRNA that generates these proteins, are crucial to the overall process. Investigating the influence of the parameters on the control target is an effective tool to reduce the model complexity and assist control system design. In order to reduce the model, we should know whether the output significantly changes under a small change of the parameters, as some parameters may not be as effective as others. By understanding how small changes in parameters affect the overall system, we can determine the most important parameters involved in generating the desired output. In the case of GAL1 mRNA, a small change in the value of the GAL4p parameter will have a large effect on GAL1 mRNA expression, indicating that this parameter significantly influences the uncertainties of GAL1 mRNA expression. Thus, this parameter is crucial to control the system. Small changes in a parameter that was not as closely linked to GAL1 mRNA expression specifically may not have a large effect on the expression. In fact, most biological systems may be sensitive to certain parameters over a large range. If the sensitivity of a parameter is small compared to other parameters, the parameter could be arbitrarily chosen from a wide range of values without altering the dynamic performance of the system. This practice has been used successfully in control systems engineering, including chemical process control, robotics control, aerospace control, etc.
For the above model, we will first conduct local sensitivity analysis that will quantify the change of the system states due to constant perturbation of a parameter at time zero. Consider a general form of nonlinear ordinary differential equations:
is the state vector with n
is the parameter vector with m
components. The initial condition of the above equations is set to x0
. Define the sensitivity coefficients of the state xi
with respect to parameter
1, …, m
and then the state's sensitivity trajectory can be described by the following:
where the state vector is denoted by
. The initial values of the sensitivity equations can be obtained by
In order to avoid scale differences in the parameters, the sensitivity coefficients are normalized as follows,
The normalized sensitivity coefficient
means that a 1 percent change of state xi
results in a 1 percent change of the parameter
over the time course t
. The sensitivity trajectory
can be calculated numerically by equations (25–26). The importance of the parameters can thus be distinguished by the magnitude of the normalized sensitivity coefficients in the time-course of the system following the constant perturbation at time zero.
The normalized sensitivity coefficients for the original model are shown in , where each small block describes a normalized sensitivity for a certain parameter at a certain time. In , the galactose input is kept at a constant,
(molec.) with the measured system output represented by m1
. It is easy to observe that the external galactose input is positively correlated to the output of m1
. The glucose network remains in the steady state when the internal glucose is emptied and there is no external glucose supply. The steady state of
13586 (molec./cell) and
3397 (molec./cell) can be determined by solving
. The function
in equation (19) was found to equal 0.9959. Thus, the glucose network can be separated from the galactose network. This is in agreement with experimental data, where a lack of glucose as a carbon source and the presence of galactose will induce the GAL
network. Based on the sensitivity analysis of m1
with respect to all the parameters (), it becomes clear which parameters can be omitted from the system. The heat map in shows the normalized sensitivity of GAL
1 mRNA concentration (m1
) after changes in all 35 parameters. The shows the sensitivities when
The represents the sensitivities when
The sensitivities of the parameters
are much smaller than those of the other parameters (near zero, green color). In general, when the external galactose
is higher, the parameters are more sensitive, with the exception of α1
Normalized parameter sensitivity trajectory with respect to time.
However, it would be careless to omit these parameters without further investigation. For example, the degradation rate
is important for achieving the steady states of the original model, so it must be retained.
come from the cooperative fractional saturation function
in the function contribute to the activation of mRNAs m1
, and m80
. However, the contributions may not be equal. If we define
, then we have
When 1 + Q >>1/P
, the above equation (28) can be approximated by
From (29), we observed that the output m1
is positively correlated to the external galactose input. The sensitivity analysis has told us that
can be set to be large enough to satisfy 1 + Q >>1/P
. Thus, only
contributes to the cooperative fractional saturation function. Similarly, the other two cooperative functions can be simplified as follows,
Therefore, the cooperative fractional saturation effect of g4d
can be approximated by the contribution from
alone. This can also be explained by the sensitivity analysis in , , and , which shows that
as well as
are less “essential” to the system output. Thus, we can set
to be zero, since they are small numbers and have small sensitivity coefficients. This results in a removal of the loops involving the GAL
, and the GAL
80p dimer complex,
. As shown in , the subsystem (
) can be simplified as
only, because the rest of the subsystems will not affect the steady state of the output. The goal of the proposed control is to maintain a level of mRNA expression at some predetermined level. The model reduction will thus have little influence on the control performance in terms of output deviation after a long enough time. The contributions from the glucose network to m1,
, are small when no glucose exists, which is necessary to induce the GAL
gene network, as explained earlier. This is consistent with the small sensitivity coefficients with respect to
. We thus set
equal to zero and
equal to a large number. Also, small sensitivity of the coefficient with respect to
leads to a slight model reduction on
Local sensitivity analysis at the steady states of four GAL mRNAs with respect to all model parameters.
Global sensitivity analysis at the steady states of four GAL mRNAs with respect to all model parameters under constant control
Thus, a Reduced Order Model can be achieved based on the analysis above. The Reduced Order Model in a state space form can thus be described as follows:
All parameters are summarized in .
The comparison of the outputs of the Original Model and the Reduced Model is shown in . The differences between the two models occur in the transition process before the first 200 min, while the difference decreases significantly in the steady state. According to the experimental results from 
, the galactose network takes 4–7 hours (240–420 minutes) to reach a steady state. Our simulation results in are consistent with the experimental data. When the external galactose level is set at
, the steady-state deviation of GAL
1 mRNA, m1
(molec./cell) is similar between the full and reduced models. However, the deviation of the transition process, the time before achieving the steady state, between the two models with the smaller external galactose level is significant. This result is consistent with the sensitivity measure in .
Comparison of original model and reduced model.
Biologically, the reduced model varies from the complex model only within the first 200 minutes of expression for relatively low galactose levels, 400 molecules. The steady-state levels of GAL1 mRNA expression after the first 200 min for both galactose levels tested were the same in both models (). In this way, the reduced model can be used to explain the expression of GAL1 mRNA, once the expression level has reached a steady state. Before the steady state is reached, only a complex model can completely model the effects of GAL1 mRNA expression from the introduction of low levels of galactose. A similar technique can be applied if the control target changes to any of the other mRNAs or proteins. This paper exemplifies GAL1 mRNA as a primary control target, but the applications of the method are not limited to a certain mRNA. Choosing GAL2 or GAL3, as well as GAL80, will produce very similar reduced models; however, using GAL4 as a control target will obviously yield a different reduced model because of the simplification in the GAL1 output system obtained by the removal of the GAL4p subsystem. A reduced model for GAL4 would have to include the GAL4p subsystem, and thus a much different reduced model would be generated.
describes the local sensitivity analysis at the steady states of four GAL
mRNAs with respect to all model parameters. The bars denote the range of the sensitivity coefficients of four GAL
mRNAs with respect to each parameter. The parameters
are smaller than 0.25, while the others are greater than one. This implies that the reaction terms associated with the insensitive parameters can be omitted no matter which mRNA is chosen as the control target. The design of a control model to control steady-state expression of GAL
3, and GAL
80 mRNA (m1
) is described in the following section.
The local sensitivity analysis describes the change of system dynamics at time t
after a step-wise perturbation to a parameter occurs at time zero. One might wonder what the effect would be when all parameters are perturbed simultaneously. Global sensitivity analysis investigates the effect of simultaneous large variations of all parameters on the states. The common methods to test global sensitivity include the Fourier amplitude sensitivity test (FAST) 
, extended FAST 
, Sobol's method 
, the partial rank correlation coefficient (PRCC) 
, as well as the weighted average local sensitivity method 
FAST is one of the classical methods to determine global sensitivity. Assuming that the time dependence of each uncertain parameter
is associated with a frequency
, we can describe that parameter by a transformation function of a sinusoid, i.e.
for all m
parameters, where s
is the time.
The first order global sensitivity of xi
with respect to the variation of the parameter
can be defined as the ratio of the variance of
to the total variance
where the total variance is
and the variance of
Although the above local sensitivity analysis pointed out some candidate parameters and their associated terms that may be removed, it only guarantees the scenario when the parameters are tuned one by one. When the parameters are simultaneously adjusted, some may become sensitive. Thus, we must also apply global sensitivity analysis to rule out corresponding parameters.
As shown in , the bars denote the range of the sensitivity coefficients of four GAL
mRNAs with respect to each parameter under constant control
The blue solid line stands for the global sensitivity coefficients for GAL
1 mRNA. In particular,
are nearly zero. The maximum of the coefficients
is smaller than 0.2, while that of the others is larger than 0.2, consistent with the local sensitivity coefficients in . This implies that the reaction terms associated with insensitive parameters can be omitted no matter which mRNA is chosen as the control target. Thus, both local and global sensitivity analyses indicate that the reduced model, equations (32–43), can be used for regulating all four GAL
As discussed earlier, any of the four GAL
mRNAs can be chosen as the control target. For demonstration purposes, we have chosen to control GAL
1 mRNA (m1
, essentially maintain the GAL
1 mRNA at a desired level
In order to keep equation (51) true, the steady state of the system becomes
As a test case for this model, we will demonstrate the control of GAL1 mRNA expression by controlling the external level of galactose. It is important to note that any gene or protein represented in the reduced model can be controlled in a similar manner. Following this same pattern, recombinant genes regulated by GAL1 could be controlled in the same manner.
Similarly, we can find the steady states for
Observing the galactose network, we can obtain two auxiliary equations:
Equation (60) indicates that the total mass of
, is up-regulated by m3
. Biologically, this is a simple concept to understand, as expression of GAL
3 mRNA, m3
, is necessary for the generation of GAL
, and the Gal3p dimer,
, in a dose-dependent manner. Similarly, the GAL
3p dimer complex is dependent on the expression of GAL
3p. Equation (61) indicates that the total mass of
, is up-regulated by m80
, in a manner similar to that described above. Both of these equations show the interplay between the biological mechanisms responsible for their creation.
Through equations (60) and (61) we can determine the remaining steady states
using the auxiliary equations
Given our desire to enforce the control condition of
, we need to control the external galactose, ,
which at steady state can be calculated as follows,
Similarly, if the control objective is to regulate GAL
The rest of the terms can be achieved using the same equations (54)–(67). Because this system has only a single steady state, one can always find a nominal control value for external galactose
to maintain any control objective at a desired level. In other words, we can always compute the necessary amount of galactose offline in order to maintain GAL
1 mRNA or any other factor, such as GAL
2 mRNA or GAL
3p, at a level of interest.
Thus, offline, we can calculate an open-loop constant control, which is the simplest type of controller that does not take into consideration state information that is fed back into the system for control purposes.
for any of the targeting mRNA, m1
, and m80
, if these mRNAs are not measurable in real-time. However, if they are measurable, one can improve the control algorithm by introducing a closed-loop feedback control. A simple proportional output feedback control for m1
can be described as
where the feedback gain
. The gain k
directly affects the convergence rate of the controlled GAL network. However, the external galactose amount should not exceed the maximum concentration of galactose. Thus, the control dosage should be bounded by
Of course, one can use a proportional-integral-derivative (PID) control in practice. Integral control can adjust the static error, but too much integral gain can reduce the stabilization of the system. Derivative control can adjust the convergence rates, however, it is too sensitive to high frequency noise. For pursuing different requirements, many more complicated control methods can be applied, such as optimal control for minimizing the galactose dosage and robust control for reducing the uncertainty, among others. Because the control objective of regulating the mRNA concentration to a desired value, is a relatively simple task, without strict requirements for reaction time and galactose dosage, a simple proportional control is adequate to accomplish the control goal.
Control Stability Analysis
Next, we will prove the stability of this equilibrium because the equilibrium is a unique steady state. If it is stable, then the state will tend to this steady state given the nominal control. The control objective m1
will be maintained at the desired level
. Biologically, this system is stable based on experimental results 
. The following is a theoretical proof based on the mathematical model.
The Jacobian matrix at this equilibrium is
, one finds the eigenvalues of the Jacobian matrix are [−895.99, −194.68, −12.54, −12.04, −0.0295, −0.0387+ 0.0024i
, −0.0387− 0.0024i
, −0.0360, −0.0015, −0.0033, −0.0033, −0.0033]. A standard stability analysis of the dynamic systems finds that the equilibrium is asymptotically stable. Thus the system will tend to the steady state asymptotically.
Because most of the eigenvalues are small in magnitude, it implies that the convergence of this controlled system will take a long time. Biologically, after we set the external galactose level to a constant value, the yeast galactose utilization network will slowly achieve balance (a steady state) and the GAL1 mRNA level will gradually reach and maintain the desired level. At this steady-state level, there will not be oscillations in the mRNA or protein levels associated with this system.
In order to change the eigenvalues to negative, we follow a common practice in control systems engineering and introduce a linear output feedback control
. The eigenvalues of the Jacobian matrix indicate how fast the control can reach a steady state. In general, the more negative, the faster. In the above control, the eigenvalues of the Jacobian matrix (73) move far away from the y axis: [−895.99, −194.67, −16.40, −5.33, −2.91, −0.026, −0.0361, −0.0360, −0.0033, −0.0033, −0.0033]. This implies that the system will converge faster under the linear feedback control. Increasing the value of k
will lead to faster convergence. However, the growth of the convergence rate is limited due to the maximum concentration of galactose. From , we can find that the internalized galactose level increases with rising external galactose levels, and the internalized galactose level enhances GAL
4, which in turn promotes GAL
1 mRNA expression. Increased GAL
1 mRNA expression will also increase the degradation rate of the internalized galactose, thus the internalized galactose will reach a steady state faster. In essence, the introduction of the feedback control loop actually increases the external galactose level, but dynamically adjusts the amount such that the GAL
1 mRNA can still reach the desired level. In this process, the degradation rate of the internalized galactose is largely enhanced, consistent with a significant increase of the magnitude of several negative eigenvalues.
We set the desired level of GAL
1 mRNA as m1
20. From equation (67), the nominal control values for [gal
* can be obtained as
molec./cell. Application of the nominal control should cause the system to slowly reach a stable equilibrium. As shown in , it takes more than 1500 min for m1
to reach the desired level. From , we can see that achieving the desired high GAL
1 mRNA level usually demands a much longer time than is required for low GAL
1 mRNA levels. The metabolic reactions involved in promoting mRNA expression and subsequent translation into proteins require a lengthy time. Multiple positive and negative feedback loops inside the network actually prevent mRNAs and the associated proteins from rapid changes due to external stimuli.
Comparison between constant and feedback control for the Reduced Order Model.
The effects of protein concentration profiles by feedback control on GAL
1 mRNA expression on the GAL
network are illustrated in . All states (unit: molec./cell) converge to a stable steady state, while g2
have slower convergence rates than the other states. It is not surprising to see that all concentrations tend to a steady state, since the system has a unique stable steady state. But the transient response of the galactose system is equally important, because a large fluctuation is not desired in a reliable bioreactor. The figure shows that the system does not have large overshoots and is thus safe for implementation in practice.
Effect of protein concentration on mRNA expression by feedback control on the GAL network.
compares the set point regulation performance of the feedback control between the original model and the reduced model using feedback control. Setting three desired levels of m1 as 10, 15, and 20, feedback control based on both the original model (solid) and the reduced model (open) can be maintained at the desired levels within a similar time frame. The reduced model-based feedback control, however, leads to slightly larger magnitudes of oscillations at the transient response. The difference is small enough not to result in any severe fluctuation in the biological system. Thus, the control design based on the reduced model can be applied to the original model directly.
Comparison of referenced tracking between the original model and the reduced model using feedback control.
The steady-state sensitivity describes the change of the steady state given a small perturbation of a certain parameter 
. It can be given by the first-order partial derivative of a state with respect to the parameter. Using the finite difference method, the normalized steady-state sensitivity (NSSS) coefficient for the steady state of state xi
with respect to the parameter
can be approximated by
is the steady state of the state xi
A small magnitude of NSSS of state xi
with respect to the parameter
implies that the parameter uncertainty of
has no significant influence on the state xi
locally. In contrast, a large NSSS implies serious influence by the parameter uncertainty. It is necessary to emphasize “locally” because the steady-state sensitivity is a local sensitivity measurement that perturbs one parameter and fixes the other parameters. As shown in , we calculated the NSSS of all states with respect to the output m1
for both constant control and the feedback control algorithm. The NSSS of the feedback control case has significantly smaller magnitude in terms of average maximum sensitivities for all states (0.93 v.s. 3.42) and output state (0.02 v.s. 1.66). Especially for the output state, the magnitude of NSSS was reduced 100-fold. This implies that the parameter uncertainty does not affect the steady state of the output. The feedback control law (81) can always drive the output to the reference level subject to uncertainties and perturbation in the parameters locally.
Comparison of normalized steady-state sensitivities (NSSS) between constant and feedback control.
Biologically, we can conclude from the above data that the feedback control is less influenced by disturbances within the system, as shown in . For mRNA expression, this means that a feedback control prevents the system from being overly sensitive to any specific parameter. The perturbations or parameter uncertainty in any parameter will not cause a significant deviation to the steady mRNA level. This is necessary for complex biological systems, where parameter uncertainty and perturbations widely exist. If the system was overly sensitive over some parameters, then it would be difficult to maintain the mRNA at a desired level. A series of checks and balances are always present to prevent the system from fluctuating out of control.
This paper employs control theory and a systems engineering approach to demonstrate the regulation in silico of the budding yeast, Saccharomyces cerevisiae, into a desired gene expression and/or metabolic activity pattern. The GAL network is used as an example to demonstrate the effectiveness of a systems approach for complex cellular system control. We demonstrate that although the structure and computational models are complicated, it is possible to use a “simple” control algorithm to manipulate the system. According to the local parameter sensitivity analysis, some terms can be eliminated from the model without adversely affecting its long-term behavior, and the control design based on the reduced model can still be applied to the original model. Both constant control and feedback control laws have been designed for the galactose network. To evaluate the control performance, we conducted local steady-state sensitivity and global sensitivity analysis. The results accordingly imply that feedback control significantly suppresses the parameter uncertainties. This sheds light on the potential to provide novel control components to eukaryotic transcriptional regulation networks, as well as other gene regulatory networks.
The approach we propose can not only be used to control the yeast GAL network, but the fundamental principles can be applied to a wide range of biological network control. Complex networks involving cell-cycle controls, DNA repair, and other genes of interest can be controlled in a similar manner, given an effective means of monitoring the output of the system. The most attractive application is the control of bioreactors. If done properly, the ability to yield a large amount of the target in a relatively short period of time, while minimizing the effects of toxic levels of gene expression products can be achieved. If these systems are tightly controlled at the cellular level, proteins of interest, such as immunogenic proteins used in the production of vaccines, restriction enzymes, biopharmaceuticals, biochemicals can be tightly regulated. This would allow for a cost-effective means of producing a large variety of products for research and commercial use. Applications to other systems include bacterial expression systems and even mammalian cell systems, based on the effectiveness of the bioreactor and the access to the requisite control variables. Of course, other physical and chemical parameters would be necessary to completely control a bioreactor-based system, including pH, temperature, and nutrient conditions. The major difference between these and the above approach is that the proposed approach focuses on cellular system control, which is still a widely open field. Our future research will implement this control system with a fully established microfluidic device to monitor the real-time expression of GAL-GFP fusion proteins. By validating the above control in an experimental environment we will demonstrate the usefulness and validity of our approach.
We are indebted to Allison Price for her editorial assistance.
1. Bhalla US, Iyengar R. Emergent properties of networks of biological signaling pathways. Science. 1999;283:381–387. [PubMed] 2. Gardner TS, Cantor CR, Collins JJ. Construction of a genetic toggle switch in Escherichia coli. Nature. 2000;403:339–342. [PubMed] 3. Cantone I, Marucci L, Iorio F, Ricci MA, Belcastro V, et al. A yeast synthetic network for in vivo assessment of reverse-engineering and modeling approaches. Cell. 2009;137:172–181. [PubMed]
4. Menolascina F, di Bernardo M, di Bernardo D. Atlanta, GA, USA: Proceedings of the 49th IEEE Conference on Decision and Control; 2010. Design and implementation of a feedback control strategy for IRMA, a novel synthetic gene regulatory network.
5. Menolascina F, di Bernardo M, di Bernardo D. Automatica; Analysis, design and implementation of a novel scheme for in-vivo control of synthetic gene regulatory networks. (In press)
6. Horak J, Wolf D. Catabolite inactivation of the galactose transporter in the yeast Saccharomyces cerevisiae: ubiquitination, endocytosis, and degradation in the vacuole. J Bacteriol. 1997;179:1541–1549. [PMC free article] [PubMed] 7. Yano K, Fukasawa T. Galactose-dependent reversible interaction of Gal3p with Gal80p in the induction pathway of Gal4p-activated genes of Saccharomyces cerevisiae. Proc Natl Acad Sci U S A. 1997;94:1721–1726. [PubMed] 8. Diep CQ, Tao X, Pilauri V, Losiewicz M, Blank TE, et al. Genetic evidence for sites of interaction between the Gal3 and Gal80 proteins of the Saccharomyces cerevisiae GAL gene switch. Genetics. 2008;178:725–736. [PubMed] 9. Johnston M, Flick JS, Pexton T. Multiple mechanisms provide rapid and stringent glucose repression of GAL gene expression in Saccharomyces cerevisiae. Mol Cell Biol. 1994;14:3834–3841. [PMC free article] [PubMed] 10. de Atauri P, Orrell D, Ramsey S, Bolouri H. Evolution of ‘design’ principles in biochemical networks. Systems Biology, IEE Proceedings. 2004;1:28–40. [PubMed] 11. Ramsey SA, Smith JJ, Orrell D, Marelli M, Petersen TW, et al. Dual feedback loops in the GAL regulon suppress cellular heterogeneity in yeast. Nat Genet. 2006;38:1082–1087. [PubMed]
12. Bennett MR, Pang WL, Ostroff NA, Baumgartner BL, Nayak S, et al. Nature; 2008. Metabolic gene regulation in a dynamically changing environment.
13. Bornholdt S. SYSTEMS BIOLOGY: Less is more in modeling large genetic networks. Science. 2005;310:449–451. [PubMed] 14. de Atauri P, Orrell D, Ramsey S, Bolouri H. Evolution of ‘design’ principles in biochemical networks. IET Sys Bio. 2004;1:28–40. [PubMed] 15. Chen KC, Calzone L, Csikasz-Nagy A, Cross FR, Novak B, et al. Integrative Analysis of Cell Cycle Control in Budding Yeast. Mol Biol Cell. 2004;15:3841–3862. [PMC free article] [PubMed] 16. Campbell RE, Tour O, Palmer AE, Steinbach PA, Baird GS, et al. A monomeric red fluorescent protein. Proceedings of the National Academy of Sciences of the United States of America. 2002;99:7877–7882. [PubMed] 17. Trumbly RJ. Glucose repression in the yeast Saccharomyces cerevisiae. Molecular Microbiology. 1992;6:15–21. [PubMed] 18. Wikswo JP, Prokop A, Baudenbacher F, Cliffel D, Csukas B, et al. Engineering challenges of bioNEMS: the integration of microfluidics, micro- and nanodevices, models and external control for systems biology. IEE Proceedings Nanobiotechnology, 2006;153:81–101. [PubMed]
19. McLean JA, Schultz JA, Woods AS. Cole RB, editor. Ion Mobility-Mass Spectrometry for biological and Structural Mass Spectrometry. Electrospray and MALDI Mass Spectrometry: Fundamentals, Instrumentation, Practicabilities, and Biological Applications: John Wiley and Sons. 2010. pp. 411–439.
20. Saltelli A, Tarantola S, Chan KPS. A quantitative model-independent method for global sensitivity analysis of model output. Technometrics. 1999;41:39–56.
21. Saltelli A, Tarantola S, Campolongo F. Sensitivity analysis as an ingredient of modeling. Statistical Science. 2000;15:377–395.
22. Sobol IM. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math Comp Simul. 2001;55:271–280.
23. Draper N, Smith H. New York: Wiley; 1981. Applied regression analysis.
24. Bentele M, Lavrik I, Ulrich M, Stosser S, Heermann DW, et al. Mathematical modeling reveals threshold mechanism in CD95-induced apoptosis. J Cell Biol. 2004;166:839–851. [PMC free article] [PubMed]
25. Varma A, Morbidelli M, Wu H. Cambridge, UK: Cambridge University Press; 1999. Parameteric sensitivity in chemical systems.