Modelling and Control of Gene Regulatory Networks for Perturbation Mitigation

Synthetic Biologists are increasingly interested in the idea of using synthetic feedback control circuits for the mitigation of perturbations to gene regulatory networks that may arise due to disease and/or environmental disturbances. Models employing Michaelis-Menten kinetics with Hill-type nonlinearities are typically used to represent the dynamics of gene regulatory networks. Here, we identify some fundamental problems with such models from the point of view of control system design, and argue that an alternative formalism, based on so-called S-System models, is more suitable. Using tools from system identification, we show how to build S-System models that capture the key dynamics of an example gene regulatory network, and design a genetic feedback controller with the objective of rejecting an external perturbation. Using a sine sweeping method, we show how the S-System model can be approximated by a linear transfer function and, based on this transfer function, we design our controller. Simulation results using the full nonlinear S-System model of the network show that the synthetic control circuit is able to mitigate the effect of external perturbations. Our study is the first to highlight the usefulness of the S-System modelling formalism for the design of synthetic control circuits for gene regulatory networks.


Introduction
In complex engineering networks such as transportation systems, power grids, irrigation networks, etc, the presence of external perturbations can have serious adverse effects on the functioning of the overall system. These undesirable effects include gridlock in the movement of vehicles, major power outages in residential and industrial areas, and unreliable water supply to farming areas. In view of this, the problem of developing a comprehensive theory of network control, particularly in the presence

Example Gene Regulatory Network
The DREAM in silico gene regulatory network challenge was established to serve as a benchmark to assess different proposed approaches to infer gene regulatory networks from given experimental data [12,13,14]. Typically, time-series data for each gene (or node) in the network are provided and the aim is to infer the underlying network, i.e. identify interconnecting edges, the direction of information flow, etc. The provided gene regulatory networks are typically subsets of actual transcriptional networks in model organisms such as E. coli and S. cerevisiae, and hence are representative of real biological systems.
In this paper, we choose the DREAM3 Size 10 data set (hereafter we use the term DREAM3 to denote this network), which consists of mRNA temporal data on a network composed of 10 interconnecting genes that is a subset of a S. cerevisiae gene regulatory network. As the dataset does not include separate protein data, in the following, we make the following two assumptions: (i) the temporal evolution of the protein is similar to the mRNA and (ii) the protein is linearly translated from mRNA. Following these two assumptions, we can lump the protein dynamics into the transcription rate of the mRNA at steady state, and this results in a complete network that can be described solely using mRNA levels. In this DREAM3 data set, information regarding the interconnectivity between each gene is provided, while the regulation type (i.e. activatory or inhibitory) is unknown. The depiction of these interactions is shown in Fig. 1(A). To facilitate the controller design procedure, a model describing the dynamics of the DREAM3 network is required, and in the following section, we discuss the selection of an appropriate modelling formalism for the DREAM3 gene regulatory network.

Michaelis-Menten and Hill-type models
Model structures employing Michaelis-Menten and Hill-type nonlinearities are commonly used to describe the dynamics of gene regulatory networks. If the regulation type and the cooperative binding are known, the modeller can either specify for an activation type of regulation or for an inhibition type of regulation. In both Eqns. (1) and (2), N P is the transcription factor, k 0 and K M are associated with the Michaelis-Menten constants and h is the Hill coefficient.
In the context of network inference, this type of model structure can be used only if the type of regulation (activatory or inhibitory) between each gene in the network is known. In the event that the type of regulation is unknown, then this model structure is not suitable as the structure of an activation or an inhibition type of regulation is different and arbitrarily assigning them in the model building stage could thus lead to poor model accuracy.  [15] whose dynamics are represented using Michaelis-Menten kinetics and Hill-type nonlinearities. For this illustration, the controller, K is a simple proportional-integral (PI) controller with the controller gains, K p = 0.01 and K I = 0.02. The pathway highlighted in yellow indicates the series of regulations involved from the control action, U to the output gene, N 3 A more fundamental problem in the context of synthetic biology is that models of this type are often not suitable for subsequent use in the design of synthetic controllers. For example, let us consider Eqn. (1) and assume that our control action (i.e. output of the controller) is given by N P . If N P K M , then F a ≈ k 0 N h P /N h P = k 0 , which renders the control action ineffective. It is thus imperative that the value of K M should be sufficiently large to ensure proper control, but as we will show below, obtaining a reliable estimate of K M from time series data is often problematic.
To illustrate the problem, we consider a model of a simple gene regulatory network taken from [15], consisting of seven interconnecting genes, as shown in Fig. 2, based on a subset of an E. coli gene regulatory network. Assume that an external perturbation enters the network through gene 1, its effect on gene 3 is measured, and fed back to a controller that regulates gene 6 through the input U. Using the standard modelling framework employing Michaelis-Menten kinetics and Hill-type nonlinearities, the associated Ordinary Differential Equations (ODEs) describing Fig. 2 are given as follows: where k 0, j , K M, j with j = 1, 2... and h are the parameters associated with the Michaelis-Menten coefficients and Hill-type nonlinearities, and γ is associated with the degradation term. Without loss of generality, for the purposes of illustration, we choose h = 1. The rest of the parameters describing Eqn. (3) are shown in Table 1. These parameters are estimated from available experimental data in [15], where one data set is used for parameter estimation and an independent data set is used for model validation. The parameters are estimated using the prediction error method with quadratic criterion, i.e.,Θ = argmin where T is the number of genes, L is the length of the data, Θ = {k 0, j , K M, j γ j } with j denotes the appropriate index describing the parameters in Eqn. (4). N i andN i represent the real experimental data and simulated data from Eqn. (3) respectively. Eqn. (4) is solved using MATLAB function fminsearch, which uses the Nelder-Mead simplex algorithm. For the controller, we choose a standard proportional-integral (PI) controller with the proportional gain, K P = 0.01 and the integral gain K I = 0.02, where these parameters can be selected using standard rules, such as the Ziegler-Nichols tuning rules (see e.g. [16]). In our simulation, shown by the solid blue line in Fig. 3, when the perturbation enters the network at time 0s it causes the expression level of N 3 to drop from its Table 1: Parameters for the network model shown in Fig. 2 using Michaelis-Menten with Hill-type nonlinearities model structure. intended reference value of 0.718 ( Fig. 3(A)). Upon sensing this drop in the expression level, the controller asserts appropriate control action, U ( Fig. 3(C)) in its attempt to bring the expression level of N 3 back to 0.718. However, as shown in Fig. 3(A), a full recovery of the output to its intended reference value is not achievable. This is because in the controller's attempt to perform the needed recovery, the exerted control action U becomes larger than K M,11 , thus the term k 0,11 U/(K M,11 +U) ≈ k 0,11 = 0.6541, which is shown in Fig. 3(D). This implies no appropriate control action can be given to the network to counter the effect of the perturbation, resulting in a large error between the output and reference value ( Fig. 3(B)). In reality, however, this may not necessarily be the case -the apparent limitation is due to the estimated value of K M,11 from the available experimental data. If the value of K M,11 is sufficiently larger than U, the 'saturation' issue is avoided. In addition, a closer look at the series of regulation along the pathway highlighted in yellow shown in Fig. 2 indicates that the values of K M, 8 and K M,13 also need to be sufficiently large in order to achieve a proper control action and recover the levels of N 3 . The problems identified above are due to the values of K M,11 , K M, 8 and K M,13 that are estimated from the available experimental data. These estimated values are relatively small when compared to the necessary control action, leading to saturated responses and large errors. Thus, a natural question arises as to whether or not these values (shown in Table 1) represent reliable estimates of these parameters. For the network shown in Fig. 2, the estimated values of K M,11 , K M,8 and K M,13 shown in Table  1 are the result of using 1 as the initial values for the parameters in the optimisation problem defined in Eqn (4). If a different set of initial values is used for the optimisation, do we obtain similar parameter values to those shown in Table 1    The plots show that the estimated parameter values are very different to the ones shown in Table 1. Using terminology from the field of system identification, there is no consistent estimate of the model parameters, as given different initial values for the optimisation, the optimiser can find different sets of parameters (see Table 2) that are equally well able to reproduce the experimental data, as shown in Figs. 4 (A), (C) and (E).

Gene Parameter Values
From Table 2, we note that there is one set of parameters that includes large values of K M,11 , K M, 8 and K M,13 . Using these larger values of K M,11 = 120.4219, K M,8 = 145.0575 and K M,13 = 99.4842, we repeat the simulation of the feedback controller shown in Fig. 2. As shown by the solid red line in Fig. 3(A), the same controller is now able to exert a proper control action to mitigate the effect of the perturbation, as the value of K M,11 is now larger than the control action, U ( Fig. 3(C)) and no issues with saturation are observed ( Fig. 3(D)).
The results shown here suggest that for this typical experimental data set and network structure, the estimated values of the model parameters, in particular K M,11 , K M, 8 and K M,13 , are not consistent. This clearly poses a significant problem when designing a controller to mitigate the effects of perturbations on this network, since different estimated values of K M,11 , K M, 8 and K M,13 lead to very different closed-loop behaviour of the control system. In light of this, coupled with the previously mentioned need for a priori knowledge of regulation type to use the Michaelis-Menten with Hill-type nonlinearities model structure, an alternative modelling formalism is clearly required in order to allow for the rational design of feedback controllers. The alternate model formalism needs to have a general structure that can accommodate both activatory and inhibitory regulations, and more importantly, the estimated model parameters from experimental show the plots using S-System model structure for genes 3, 6 and 7 respectively. The notation p 0 denotes the parameter set obtained when initial value of 1 is used for the optimisation (shown in Table 1). The notation p : 0.1, 1, 10 → p 0 indicates the estimated parameters using initial values of 0.1, 1 and 10 are similar to p 0 . data should be consistent, so that it can be reliably used for controller design.

S-System models
The so-called S-System modelling formalism has been proposed as an alternative approach to describe the dynamics of gene regulatory networks. The S-System modelling framework was originally developed from the field of biochemical system theory (see e.g. [17,18]), and when it has been used to describe the dynamics of gene regulation (see e.g. [19,20]), it has been shown to be as accurate as Michaelis-Menten with Hill-type nonlinearity models (see [21]). In particular, the authors in [21] rigorously analysed the 'validity' range of the concentrations produced by both S-System and Michaelis-Menten models to determine which model differs most from the 'true' concentration obtained via experiment. It was found that, not only were S-System models  as accurate as Michaelis-Menten type models within the same concentration range, but the S-System models were more accurate over a wider range of concentrations. Based on this and other analyses, the authors suggested that the S-System model formalism better represents the actual biochemical system. The S-System models we consider in this work have the following form: where i denotes the number of biochemical component, a i > 0, b i < 0 and c i, j ∈ (−∞, +∞) are constants, N i represents the biochemical component, M 1 and M 2 are the total number of components involved in the interaction, U j is the external input and M 3 is the number of input. The power exponent terms, p i, j and q i, j are associated with the production and degradation terms respectively. For simplicity, we assume q i, j = 1 throughout this paper, so that a positive value for the parameter p i, j represents activation while a negative value represents inhibition. Thus, the S-System model has a general structure that can accommodate either an activation or inhibition type of regulation via the sign of p i, j , and no prior knowledge of the type of regulation at each node in the network is required in the model building process. The S-System model describing the gene regulatory network shown in Fig. 2 is given as follows: Note that for dN 1 /dt, a constant value denoted by d 1 is added to the model to ensure the overall mRNA level stays positive since D is negatively correlated with N 1 and b 1 is negative due to the degradation term. Like before, we use one set of experimental data for parameter estimation and an independent set of data for model validation. The parameters are estimated using the prediction error method with quadratic criterion (Eqn. (4)) with Θ = {a i , b i , c i , d 1 , p i, j } where i and j denote the appropriate indices in Eqn. (6). The estimated parameters, using 1 as the initial value for all parameters in the optimisation, are given in Table 3.
We repeat the feedback control design using the same configuration shown in Fig.  2. The feedback control response when a perturbation enters the gene regulatory network is shown in Fig. 5. When the S-System model is used, the controller is able to produce an appropriate control action to attenuate the effect of the disturbance. There is no saturation issue observed, unlike in the scenario where the Michaelis-Menten with Hill-type nonlinearities model structure is used.
We proceed further to check whether the estimated parameters for the S-System model are consistent or not. As before, we choose the initial parameter values for the optimisation to be 0.01, 0.1, 10 and 100. The resulting estimated parameters are given in Table 4. The results shown in Figs. 4 (B), (D) and (F) indicate that, using this model structure, the estimated parameters are now consistent. Denoting p 0 as the estimated parameter set obtained when 1 is used as the initial value for optimisation, we observe that when initial values of 0.1 and 10 are used, the estimated parameters are Table 3: Parameters for the network model shown in Fig. 2 using S-System model structure.
Gene Parameter Values close to p 0 (see Table 4). When initial values of 0.01 and 100 are used, the estimated parameters are not close to p 0 , but in this case the model responses do not reproduce the experimental data. Taken altogether, these results suggest that we are able to obtain consistent estimates of the model parameters from experimental data when using the S-System model structure, making this modelling formalism much more suitable for use in the design of feedback controllers for perturbation mitigation.

Identification of an DREAM3 network using S-System model
In the previous section, we have illustrated why the S-System model formalism is a more appropriate way to model gene regulatory networks for the purposes of control system design. We now proceed to use the S-System model structure to identify, model, and design a biologically implementable perturbation mitigation controller for the DREAM3 network. Fig. 1(A) shows the interconnection between the genes in the DREAM3 gene regulatory network. In contrast to the network shown in Fig. 2, here no information is provided regarding the type of regulation between the interconnecting genes, and therefore we use system identification techniques (see e.g. [22]) to infer the type of regulation within the network. Note that, since no information regarding the type of regulation between the interconnecting genes is available, the Michaelis-Menten with Hill-type nonlinearities model structure cannot be used in this case. System identification techniques have been used to build models of gene regulatory networks in several previous studies, including [23,24,25], where linear black box network models were considered and the directions and the types of regulation were identified based on available data on gene expression profiles. In this paper, we consider a nonlinear grey box S-System model, given that we have prior knowledge about the network interconnections, and focus our attention on the identification of the type of regulation between the interconnecting genes. We use one data set for parameter estimation and another independent data set for model validation. Note that both the estimation and validation data sets used are the provided temporal profiles from the DREAM3 gene regulatory network challenge. The S-System model for the DREAM3 gene regulatory network following Fig. 1(A) is given by dN 10 dt = a 10 N p 10,1 7 Again note that for dN 9 /dt, a constant value denoted by d 9 is added to the model to ensure that the overall mRNA level stays positive since U 3 is negatively correlated with N 9 and b 9 is negative due to the degradation term. The parameters are estimated using Eqn. (4) with Θ = {a i , b i , c i , d 9 , p i, j } and T = 10.
Using 1 as the initial value for all parameter in the optimisation, the estimated parameters of Eqn. (7) are given in Table 5. Fig. 6 shows the comparison between the S-System model and the real data on the validation data set. The initial conditions for solving the ODEs are the first data points of each gene taken from the experimental data set.
From the estimated parameters shown in Table 5, we are able to determine the type of regulation in the network, where a positive value of the power term denotes activation while a negative value of the power term denotes inhibition. Reassuringly, all the a priori known degradation terms were identified to have negative values, in accordance with current biological data on the network.
The comparison between the S-System model and the real data on the validation data set shows good agreement, suggesting a good level of accuracy of the model. To  Exp. Data Model Figure 6: Comparison between S-System model and DREAM3 data on the validation data set that is not used for parameter estimation. quantify this, we calculate the Mean Square Error (MSE) for each gene between the S-System model and the real data. The MSE is computed using, where L is the length of the data, N i andN i respectively represent the experimental and the simulated data and i = 1, 2, . . . , 10. Table 6 shows the computed MSE for both the estimation and validation data sets. The total MSE, MSE T , is obtained by summing all the individual MSE from each genes. In general, the MSE values are small and similar between the two data sets. With the regulation types in the DREAM3 network as identified, the network interactions are as shown in Fig. 1(B).

Modelling of DREAM3 with Michaelis-Menten with Hill-type nonlinearties
Now that the regulation types between each node (activation or inhibition) have been identified, we can also use Michaelis-Menten with Hill-type nonlinearities to model the DREAM3 network, as follows: We want to investigate whether the Michaelis-Menten with Hill-type nonlinearities model would encounter the same problem of inconsistent parameter estimates as highlighted in Section 3.1. For the purposes of illustration, we focus only on the highlighted pathway that involves the series of regulation from the control action to the output gene (see Fig. 1(C)) and as before set h = 1.
We repeat the parameter estimation exercise (i.e., using Eqn. (4)) where we choose 0.01, 0.1, 10 and 100 as the initial values for the optimisation for both Michaelis-Menten with Hill-type nonlinearities and S-System model structures, focusing only on genes 1, 4 and 9. The results are shown in Fig. 7 and the estimated model parameters are given in Tables 7 and 8 As shown in Fig. 7, the estimated parameters using the Michaelis-Menten with Hill-type nonlinearities model are not consistent, as different sets of parameter are able to reproduce the dynamics of the experimental data equally well. For the S-System model, however, we obtain consistent estimates of the model parameters for genes 1 and 9 when the initial values used for optimisation are 0.01, 0.1 and 1, while for initial values of 10 and 100, the resulting parameters cannot reproduce the experimental data. For gene 4, we obtain consistent estimates of the model parameters when the initial values used for optimisation are 1, 10 and 100, while for initial values of 0.01 and 0.1 there is again poor agreement between model responses and experimental data.

Discussion on the parameter estimates of the model structures
Through our analysis of different modelling formalisms for the gene regulatory networks considered here, we have illustrated the inconsistent estimates of the model parameters obtained when using Michaelis-Menten with Hill-type nonlinearities model. This means that these model parameters are not identifiable from the available experimental data. One reason for this could be that these experimental data do not excite the relevant dynamics (in particular the saturation region) thus making the data not informative enough to obtain a consistent estimate. This inconsistent estimate is related to the notion of 'practical parameter identifiability' (see e.g. [26], [27]) where the available experimental data is unable to excite the relevant dynamics to provide consistent estimate for a given model structure, as observed here. The problem of inconsistent parameter estimates is also observed in [28], where the authors attempt to build a comprehensive network model for the plant circadian system, and the interactions between genes are modelled using the Michaelis-Menten with Hill-type nonlinearities model structure. The model parameters are estimated from experimental data, which are the temporal profiles of the circadian genes and proteins and a total of eight different parameter sets are found to be able to reproduce the experimental data. The estimated  Although its relevance from the point of view of control system design has not todate been considered, the problem of obtaining consistent estimates of parameters in the Michaelis-Menten model structure has been previously investigated (see the review paper [29] and references therein). In [30] and [31], different methods for fitting the Michaelis-Menten equation were analysed, and both studies concluded that different fitting methods will give different estimates of the parameters unless the experimental data is free from error (which in biological reality it never is). Different approaches to estimate the Michaelis-Menten coefficients have also been studied in [32], [33] and [34], and those studies concluded that it is difficult to obtain a consistent estimate of the Michaelis-Menten coefficients unless particular design considerations are taken into account.
On the other hand, for the parameters of the S-System model, our two illustrative examples indicate that these parameters are locally identifiable [35], as we are able to obtain consistent parameter estimate when different initial values are used for the optimisation. The identifiability of model parameters using a power law type of model structure (that includes the S-System model) has been investigated in [36]. Their analyses show that while in general it is practically challenging to obtain consistent estimate for all the parameters in the model, one can obtain consistent estimates of the model parameters under certain conditions. Recent work by [37] also shows that with an ap- propriate choice of optimiser, one can obtain consistent parameter estimates using the S-System model structure.

Design of a Feedback Controller for Perturbation Mitigation
Here, we show how the S-System model of the considered gene regulatory network can be used to design a controller for perturbation mitigation. To achieve an implementable design, a genetic-based controller is required, and there are frameworks available for such designs (see e.g. [38,39]). In this paper, we employ a frequency domain control design methodology, motivated by the design framework proposed in [39]. In order to design controllers in the frequency domain, a linear model is required. As the S-System is a nonlinear model, we linearise it to obtain a transfer function model using the sine sweeping method (see e.g. [22,40]).

Sine sweeping method
In the sine sweeping method, sinusoidal input signals over the frequency range of interest are given as the inputs to the system. The output responses within the frequency range are then analysed in terms of their magnitude and phase relative to the input signal. By collecting these magnitude and phase values, the frequency response and transfer function model of the system can be easily obtained. Here, we summarise the procedure for obtaining a transfer function model using the sine sweeping method method and refer readers to [22,40] for complete details. Consider a sinusoidal input u(t) = A sin(ω 0 t), where A is the amplitude and ω 0 is the frequency. For any linear time invariant system, the output would be also sinusoidal with the same frequency but with scaled amplitude and a phase shift. In practice, the output response is subject to transient effects, as well as the effects of nonlinearities and disturbances d(t), yielding, where B = A|G( jω 0 )|, φ = ∠G( jω 0 ) = tan −1 Im|G( jω 0 )| Re|G( jω 0 )| and G( jω 0 ) is the transfer function relating the input and output with j denotes the imaginary number.
The effects of transients and nonlinearities can be removed by neglecting the initial part of the data and assuming that the linear dynamics make the dominant contribution to the overall response. To reduce the effect of d(t) on y(t), one can use a correlation method [22], where the idea is to correlate y with a sine and cosine of the same frequency and average it over the length of the data N L (see Fig. 8).
Substituting Eqn. (10) into (11), and after some algebraic manipulation, we arrive at From Eqn. (12), the second term for both I S (N L ) and I C (N L ) will go to zero as N L → ∞. Assuming d(t) is a stationary stochastic process with zero mean value and covariance function R d (l) such that ∑ ∞ l=0 l|R d (l)| < ∞, the third term for both I S (N L ) and I C (N L ) will be zero as N L → ∞, since the variance of the third term decays at a rate of 1/N L (see [22] for details). From the remaining terms of Eqn. (12), the magnitude, |G( jω 0 )| and the phase, ∠G( jω 0 ) can be estimated using the following equations, i.e.
For the DREAM3 network, we assume that the input to the network is through U 3 and the output of interest is the expression of gene N 1 . We apply sinusoidal signals of the form 3 sin(ωt) + 3 with the frequency ω ranging from 0.001 rad/s to 1.000 rad/s. Despite using a nonlinear model, we note that the output sinusoidal responses have the same frequency as the input and no subharmonics are apparent, indicating a dominant linearity of the model. By computing the magnitude and phase values using Eqn. (13), the Bode plot of the DREAM3 network from input U 3 to output N 1 is obtained and shown in Fig. 9.
From the Bode plot, we note the following: (i) At low frequency, the magnitude of the system is about -22.5dB. (ii) The corner frequency is 0.11 rad/s. (iii) At the  corner frequency, the slope is close to -40dB/dec and the phase is approximately -90 • , suggesting a second order system with repeating poles. Thus, the transfer function relating input U 3 to output N 1 can be approximated by 0.0009 s 2 + 0.22s + 0.012 (14) From the sine sweeping method, the linear transfer function of the gene regulatory from U 3 to N 1 is given by Eqn. (14). We compare the accuracy of the linear model with the nonlinear S-System model through a step response comparison, as shown in Fig. 10. Since the base signal level used in the sine sweeping method is 3, the input is stepped from 3 to 4.
From Fig. 10, we observe similar performance between the two models in terms of their transient responses, i.e. similar rise time and settling time. On the other hand, the steady state levels between the two models are different with the linear model having a higher steady state level compared to the nonlinear model. Nevertheless, the difference between these two steady state level is relatively small, indicating acceptable accuracy of the linear model in approximating the nonlinear S-System model relating input U 3 to output N 1 .
With this transfer function identified, we can proceed with the design of the controller using a frequency domain approach.

Design of a genetic phase lag controller
Here, we illustrate the design of the genetic phase lag controller. A phase lag controller is chosen, as this type of controller is typically used to improve disturbance rejection and reduce steady state errors. The phase lag controller has the following form: where the zero of the controller z = −(a P + (K 1 /K 2 )) and the pole of the controller p = −a P , with the gain of the controller being K 2 . As both the gain and phase margins of the system obtained from the Bode plot are infinite, our primary focus is on improving the transient dynamics of the disturbance rejection and reducing the steady state error. The transfer function given in Eqn. (14) is a type 0 system, and with the use of a phase lag controller, there is no integrator in the open loop gain to eliminate the steady state error. As such, when choosing the pole of the phase lag controller, we try to place the pole, a P as close as possible to the origin. Likewise, the static error constant, K p = 0.0027K 2 should be chosen as large as possible to reduce the steady state error. The choice of the design parameters are constrained by the achievable biological values and following the range of allowable values given in [39]; the following allowable parameter ranges are adhered to: 0.0002 ≤ a P ≤ 0.0040, K 1 < 2.3 and K 2 < 1.8.

Simulation Results
While the design of the controller is carried out using the linear model, for implementation, we carried out our simulation using the nonlinear S-System network model. In most gene regulatory network perturbation mitigation problems, we are interested in maintaining the steady state level of a particular gene of interest in the presence of a perturbation. Biologically, this can be interpreted as maintaining the level of expression of a gene of interest to ensure optimal biological function. Thus, in this simulation example, we are interested in maintaining the steady state level of N 1 at its desired reference value in the presence of a perturbation. Here, we assume that the perturbation enters the network through U 1 and our control action is provided by U 3 as depicted in Fig. 1(C).
In the absence of a perturbation, the steady state level of N 1 is 0.486, thus, our control objective is to maintain the steady state level of N 1 close to 0.486 in the presence of a perturbation. In our simulation, a perturbation in the form of a step response with amplitude of 2 enters the network at time 4000s. As can be seen in Fig. 11(A), without control, the steady state level of N 1 increase to 0.63 and is unable to return to its desired value.  In the design of the phase lag controller, the following values are chosen. To have the pole close to the origin, we choose a P = 0.0002. To have the static error constant as large as possible, we choose K 2 = 1.7. For K 1 , we initially consider two cases, i.e. K 1 = 0.04 (controller's zero close to origin) and K 1 = 2 (controller's zero far from the origin). The simulation results are shown in Fig. 11(B). For a small value of K 1 , we see that the performance of the system is slow and at time 6000s, there is still a noticeable steady state error, i.e. 0.044. On the other hand, for a large value of K 1 , we see a significant improvement in the performance, where we get a faster response and an almost zero steady state error, i.e. 0.0008.
The Bode plots of the system with and without control are shown in Fig. 12. For a small value of K 1 , we note that the phase margin of the system is 97 • . On the other hand, for a large value of K 1 , despite the good performance, we note that the phase margin of the system reduces from 97 • to 10 • , which is less than typically specified values. Thus, a compromise between the transient performance and overall stability robustness needs to be performed when designing the controller, and this trade-off can be effectively managed through the choice of the controller parameter K 1 . According to standard specifications, the phase margin is typically required to be between 45 • to 60 • (see e.g. [16]) to achieve satisfactory performance. To find the 'optimal' value of K 1 that can achieve fast response, small steady state error and achieve a phase margin in the aforementioned range, we proceed as follows.
The transfer functions of the process and the lag compensator are given by Eqns. (14) and (15) respectively. Rewriting them here together with the substitution of a P = 0.0002 and K 2 = 1.7, as well as defining G OL (s) as the open loop gain transfer function, we have the following expression.
As shown by the green dash-dotted line in Fig. 11(B), with K 1 = 0.8, the magnitude plot has shifted to the left. This left shift in magnitude changes the gain cross over frequency from 0.1 rad/s to the one we specified, i.e. 0.05 rad/s. On the other hand, the phase plot is similar to the case when using large K 1 . Nevertheless, more importantly, the Bode plot shown in Fig. 12(A) and (B) shows that the new phase margin is 47.4 • when using K 1 = 0.8, which is within the preferred range and a significant improvement compared to using large K 1 .

Conclusions
Although several modelling formalisms are now available for the representation of gene regulatory networks, the question of their suitability for the design of synthetic feedback control systems has so far received little attention in the literature. In this paper, we show that standard modelling approaches employing Michaelis-Menten models with Hill-type nonlinearities are not appropriate for use in the design of synthetic controllers, for two reasons. Firstly, such models require the type of regulation between interacting genes in the network to be known a priori, which is highly unlikely to be the case in general. Even more problematically, the values of the particular parameters in such models on which the controller design depends cannot in general be reliably identified from standard time-series data.
As an alternative approach, we propose the use of the S-System modelling formalism. While the use of the S-System modelling formalism for describing the dynamics of gene regulatory networks is well established, its usefulness for the purposes of control design has not so far been investigated. Here, we showed that using this modelling formalism combined with standard system identification procedures allows us to establish the type of regulation between each gene, obtain consistent estimates of model parameters, and hence derive a model that is suitable for the design of a synthetic genetic feedback controller. Given that the design of the considered genetic feedback controller is carried out in frequency domain, we showed that the nonlinear S-System model can be approximated by a second order linear transfer function using the sine sweeping method. Based on this transfer function model, we designed a genetic phase lag feedback controller, whose structure and parameter values can be readily implemented biologically. Simulation results show satisfactory performance of the controller in mitigating external network perturbations. The proposed modelling and control system design approach considered here has been tailored to the problem of mitigating external perturbations in gene regulatory network. However, the proposed approach can be readily extended to address other control problems (e.g. reference tracking) and should have wide potential application to network control problems throughout the field of synthetic biology.