# Robust Genetic Circuit Design: A MixedH∞ and IQC Analysis

Authors: Ruijuan Chen; Wei Pan; Dongfei Fu; Xin He
Journal: Advances in Mathematical Physics
 (2017)

---

## Abstract

This paper considers the problem of designing a genetic circuit which is robust to noise effect. To achieve this goal, a mixedH∞ and Integral Quadratic Constraints (IQC) approach is proposed. In order to minimize the effects of external noise on the genetic regulatory network in terms of H∞ norm, a design procedure of Hill coefficients in the promoters is presented. The IQC approach is introduced to analyze and guarantee the stability of the designed circuit.


---

## Body

## 1. Introduction

Genetic regulatory network (GRN) is subjected to noise disturbances that may occur at transcription, translation, transport, chromatin remodeling, and pathway specific regulation. The GRN diagrams that resemble complex electrical circuits are generated by the connectivity of mRNAs and proteins [1]. Mathematical and computational tools have been utilized to develop the genetic circuits and systems using biotechnological design principles of synthetic GRN, which involves new kinds of integrated circuits like neurochips inspired by the biological neural networks [2]. This method leads to a large-scale system composed of several interconnected subsystems. The previous work [3] performs a hierarchical analysis by propagating the IQC characterization of each uncertain subsystem through their interaction channels. More specifically, both plant states and the IQC dynamic states are used as feedback information in the closed-loop system model, and then the robust l2 stability analysis is performed via dynamic IQCs. Thereby, the synthesis conditions for the proposed full-information feedback controller are derived for the linear matrix inequality (LMI) systems [4].Therefore, stability analysis of uncertain GRN is a prerequisite for any design issue. From the perspectives of control engineering,H∞ is a key performance index to evaluate the noise rejection/attenuation capability. Unlike the external control inputs used in the conventional robust control theory [5], the feedback regulation mechanism is embedded in the GRN. We construct a genetic circuit by introducing a Hill function type feedback loop from proteins (mostly from transcription process) to regulate the expression of target genes. By binding to promoter domain, the GRN is mean square asymptotically stable with a given noise attenuation level γ.The paper is organized as follows. Section2 introduces the mathematical model of GRN. A design procedure for the Hill coefficients is proposed in Section 3. We provide an example to illustrate the developed design method in Section 4. The concluding remark is given in Section 5.


## 2. Problem Formulation

The activities of a gene are regulated by other genes through their interactions, that is, the transcription and translation factors [6, 7]. The underlying dynamics can be modeled as a gene i=1,2,…,n,(1)dmitdt=-limit+∑j=1nGijhjpjt+eit,dpitdt=-cipit+dimit,where mi(t),pi(t)∈R are concentrations of mRNA and protein of the ith gene at time t, respectively, li,ci∈R+ are the degradation rates of the mRNA and protein, di∈R+ is the translation rate, ei(t) is the external noise, and (2)hjx=βjx/kjnj1+x/kjnj,j=1,…,n,is a monotonically increasing function [8] in which nj is the Hill coefficient, βj is a positive constant, and kj is the apparent dissociation constant derived from the law of mass action, which equals the ratio of the dissociation rate of the ligand-receptor complex to its association rate. The family of positive Hill functions is shown in Figure 1. In this paper, the Hill function assumes that protein j is an activator of gene i [6, 7]. The matrix G=Gij∈Rn×n is the coupling matrix of the GRN and ei is defined as a base rate. System (1) can be written into the compact matrix form:(3)dmtdt=Lmt+Ghpt+et,dptdt=Cpt+Dmt,where m(t)=m1(t),…,mn(t)T, p(t)=p1(t),…,pn(t)T, L=diag-l1,…,-ln, C=diag{-c1,…,-cn}, D=diagd1,…,dn, e(t)=[e1(t),…,en(t)]T, and h(p(t))=[h1(p1(t)),…,hn(pn(t))]T. To simplify our exposition, we use a more general set of notations and shift the equilibrium point of the noiseless system to P; then model (3) can be expressed as(4)dxtdt=Axt+BHxt+Et,where(5)A=L0DC,B=0G∗0,with ∗ being an arbitrary matrix such that |B|≠0, (6)xt=mtpt,Hxt=0hpt,Et=et0.In this model, the system states of mRNAs and proteins play different roles in regulation, for example, activators, repressers, or other factors. We name x(·)=[x1(·),…,xn(·),xn+1(·),…,x2n(·)]T∈R2n as the deviation of concentration from the equilibrium point of (3). The rate of change in xi denoted by x˙i represents the concentration changes of the variables due to production or degradation. H(·)=0,…,0,h1(·),…,hn(·)T represents the regulation function on the ith variable, which is generally a nonlinear or linear function on the variables [x1(·),…,xn(·),xn+1(·),…,x2n(·)]T, but has a form of monotonicity with each variable. The degradation parameters matrix A has zero elements on its nondiagonal plane; the matrix B defines the coupling topology, direction, and the transcriptional rate of the GRN. When the input is u(t)≜Hxt, the system model (4) can be rewritten as(7)x˙t=Axt+But+Et,ut=Hxt,where E(t) is the vector of zero-mean white Gaussian noise.Figure 1
Positive Hill function.In this paper, we aim to address the following problem.Problem 1.
Given the system represented by model (7) and parameters A, B, the purpose of robust genetic circuit design is to determine the parameters in H(·) such that (i)
the whole system is stable;(ii)
H∞ norm of the noise e(t) in the measurement channel x(t) is minimized (we assume full observation).In the following section, we propose a mixedH∞ and IQC approach to tackle Problem 1. The objective of this approach is to promote H∞ method in the design of Hill function for GRN. We will give the theoretical analysis to underpin this technique.


## 3. Analysis and Design

We take pointP as the equilibrium position due to the special shape of Hill functions in Figure 1. The intuitive way to address such an issue would be firstly to design a static feedback u=Foptx by which the closed-loop system can achieve the minimum γmin. Then we prove the stability of the system with the nonlinearity involved by analysis methods [9]. However, the optimal performance with u=Foptx for the linear system does not necessarily guarantee the optimal performance for the nonlinear system with u=H(x). In some cases, the nonlinearity of u=H(x) might even worsen the system performance to a degree which is far from optimal. Therefore, robust performance has rarely been considered for the nonlinear system. In the following subsections, we resolve this difficulty.


### 3.1. Preliminaries

We first recall some preliminary results in the system analysis via IQCs from [10]. Let RH∞m×n be the set of real proper rational function matrices without right-half plane poles and let L2el[0,∞) be the set of functions f:[0,∞)→Rl that have finite energy on the interval [0,T], ∀T>0; that is,(8)f22=∫0Tft2dt<∞,∀T>0.The Fourier transform for an Rl-valued function f:[0,∞)→Rl is denoted as f^(jω). Consider the feedback configuration in Figure 2,(9)v=Gw+f,w=Δv+e,where f∈L2em[0,∞), e∈L2en[0,∞), and G and Δ are two causal operators. Note that G is stable and Δ is bounded but could be nonlinear, time-varying, or uncertain. Δ is said to satisfy the IQC defined by Π if the two vectors of signal v, w fulfill(10)∫-∞+∞v^jωw^jω∗Πjωv^jωw^jωdω≥0,in which Π:jR→C(m+n)×(m+n) can be any measurable Hermitian-valued function defined on the imaginary axis and the superscript ∗ denotes the complex conjugate transpose.Figure 2
The feedback system model structure.Lemma 2 (see [10]).
LetG(s)∈RH∞m×n and Δ be a bounded causal operator. If the following assumptions hold, (i)
for everyτ∈[0,1], the interconnection of G and τΔ is well-posed;(ii)
for everyτ∈[0,1], the IQC defined by Π is satisfied by τΔ;(iii)
there existsϵ>0 such that(11)GjωI∗ΠjωGjωI≤-ϵI,∀ω;then the feedback interconnection of G and Δ is stable.Remark 3.
From Lemma2, an important conclusion on IQC stability analysis can be drawn: if Δ satisfies several IQCs, it also satisfies all nonnegative linear combinations of IQCs. For example, if v, w satisfy the IQC defined by Π1 and the IQC defined by Π2, they will also satisfy the IQC defined by λ1Π1+λ2Π2 for all λ1,λ2≥0. In other words, the set of IQCs that Δ satisfies forms a “description” for the block Δ, and the more IQCs we know that Δ satisfies, the more precisely that we can describe the uncertainty of Δ.Lemma 4 (KYP lemma, [11]).
SupposeM(s)=Ccl(sI-Acl)-1Bcl+Dcl. Assume (Acl,Bcl) is stabilizable and Acl has no eigenvalues on the imaginary axis. Then the following conditions are equivalent: (i)
The systemM is stable, and M∞<γ.(ii)
We have∀ω∈[0,∞),(12)jω-Acl-1BclI∗CclTCclCclTDclDclTCclDclTDcl-γIjω-Acl-1BclI<0;(iii)
There exists a symmetric matrixP>0, such that(13)AclTP+PAclPBclBclTP-γI+CclTDclTCclDcl<0.Lemma 5 (see [11]).
By using linear fraction transformation, we convert the system configuration from Figure2 to Figure 3. Let e∈L2en[0,∞), Gl(s)∈RH∞ with corresponding dimension, and let the system L2-gain be the performance measurement.
Assume thatΔ satisfies an IQC defined by Π having the following block structure:(14)Π=Π11Π12Π12∗Π22;then the system in Figure 3 is stable and has robust L2-gain γ, if (i)
for everyτ∈[0,1], the interconnection of Gl and τΔ is well-posed;(ii)
for everyτ∈[0,1], the IQC defined by Π is satisfied by τΔ;(iii)
the frequency domain inequality,(15)GljωI∗ΓGljωI<0,holds for allω∈[0,∞), where (16)Γ=I0000Π11jω0Π12jω00-γ2I00Π12∗jω0Π22jω.Figure 3
The transferred system model structure.


### 3.2. Design Procedure

The algorithm for the design of Hill function considering both the stability and performance of the nonlinear system will be proposed.Proposition 6.
After shifting the equilibrium pointP (Figure 1) to the origin, Hill function h(x) takes the following form:(17)hx=βx+knkn+x+kn-β2,and there exist real scalars α1≥α2≥0, such that h(x) satisfies the IQCs defined by(18)π1=-2α1α2α1+α2α1+α2-2,π2=β22zy-y-z,where z≥0.Proof.
π
1 is from a sector bound condition α2v2≤h(v)v≤α1v2; to see this we notice that in the time domain(19)vhvTπ1vhv=2α1v-hvhv-α2v≥0.The derivation of π2 involves an approximation, in which h(v(t))=δ(v(t))v(t) with a scalar function δ∈L∞ with δ∞≤β/2; then(20)vhvTπ2vhv=vTβ22z-δ2z+δy-yv≥0.This ends the proof.Remark 7.
From Proposition6, we should know that the set of IQCs is a sufficient condition for stability; therefore, it is conservative. Yet the challenge for nonlinear system stability is significant, because all the existing methods, including Lyapunov theory and IQCs, are conservative.In Proposition6, π1 and π2 are for the case of scalar x and h(x). Here we consider the case of vectors x and h(x).Corollary 8.
Whenx=[x1,…,xi,…,xn,xn+1,…,x2n]T and H(x)=[0,…,0,h1(x1),…,hi(xi),…,hn(xn)]T, then IQCs for H(x) take the following forms:(21)Π1=Π1–11Π1–12Π1–21Π1–22,where (22)Π1–11=diag-2α11α21,…,-2α1iα2i,…,-2α1nα2n,Π1–12=diagα11+α21,…,α1i+α2i,…,α1n+α2n,Π1–21=Π1–12T,Π1–22=diag-2,…,-2,…,-2,Π2=Π2–11Π2–12Π2–21Π2–22,where (23)Π2–11=diagz1,…,zi,…,zn,Π2–12=diagy1,…,yi,…,yn,Π2–21=-Π2–12T,Π2–22=diag-z1,…,-zi,…,-zn,with zi≥0. Combining Π1 and Π2, we obtain the following IQCs for H(x):(24)Π=Λ1Π1–11Λ1Π1–12Λ1Π1–21Λ1Π1–22+Λ2Π2–11Λ2Π2–12Λ2Π2–21Λ2Π2–22,for any n×n nonnegative diagonal matrices Λ1 and Λ2.Example 9.
Taken=1, n1=3, k=1, and β=1 as an example; therefore,(25)hx=x+131+x+13-0.5.Assume that the evolution of x is within a limited a range [-1,1]; then a sector condition of α1=0.84, α2=0.35 with δ∈[-0.5,0.5] can form a good bound on h(x), which is shown in Figure 4.Figure 4
Hill functionh(x) and its bounds. In this picture, the solid line is the considered Hill function and dashed lines are the corresponding bounds.Proposition 10.
Consider the system configuration in Figure5, let e∈L2em[0,∞), Gl(s)∈RH∞ be the genetic network, and x(·)=[x1(·),…,xn(·),xn+1(·),…,x2n(·)]T and H(x)=[0,…,0,h1(x1),…,hi(xi),…,hn(xn)]T represent the designed feedback Hill function. We denote with ΔH the uncertainty derived from H(x), which satisfies the IQCs in (24).
Then the system in Figure5 is stable and has robust L2-gain γ, if (i)
for everyτ∈[0,1], the interconnection of Gl and τΔH is well-posed;(ii)
for everyτ∈[0,1], the IQC defined by Π is satisfied by τΔH;(iii)
the frequency domain inequality,(26)GljωI∗ΓGljωI<0,holds for allω∈[0,∞], where (27)Γ=I0000Π11jω0Π12jω00-γ2I00Π12∗jω0Π22jω,andΠijs take the form in (24).Figure 5
The system with Hill feedback model structure.Proof.
The result follows from applying the IQCs defined by (24) in Corollary 8.Proposition 11.
Let the Hill coefficientsHc be the set of proper values of (β,k) for a Hill function; given any (β,k)∈Hc and n, a Hill function can be uniquely defined and hence a set of IQCs in the form of (24) can be formulated. Then the design of Hill function feedback can be reformulated into the following linear matrix inequality (LMI) problem: (28)infΛ1,2>0,Hc,ωγ,such that (26) is satisfied.The frequency-dependent inequality in (26) can be resolved by using YALMIP (https://users.isy.liu.se/johanl/yalmip/). Since the linear matrix inequality (28) has to be solved simultaneously for all frequencies in the frequency domain, the basic strategy is to start with a small set of points and then check whether the solution satisfies all ω. If it does not, add one or more frequency points to the previous set and solve LMI again. Another more efficient way to solve the LMI is to transform (26) into the time domain by using the KYP lemma of Lemma 4 if the state-space realizations of Gl and Π have been obtained.


## 3.1. Preliminaries

We first recall some preliminary results in the system analysis via IQCs from [10]. Let RH∞m×n be the set of real proper rational function matrices without right-half plane poles and let L2el[0,∞) be the set of functions f:[0,∞)→Rl that have finite energy on the interval [0,T], ∀T>0; that is,(8)f22=∫0Tft2dt<∞,∀T>0.The Fourier transform for an Rl-valued function f:[0,∞)→Rl is denoted as f^(jω). Consider the feedback configuration in Figure 2,(9)v=Gw+f,w=Δv+e,where f∈L2em[0,∞), e∈L2en[0,∞), and G and Δ are two causal operators. Note that G is stable and Δ is bounded but could be nonlinear, time-varying, or uncertain. Δ is said to satisfy the IQC defined by Π if the two vectors of signal v, w fulfill(10)∫-∞+∞v^jωw^jω∗Πjωv^jωw^jωdω≥0,in which Π:jR→C(m+n)×(m+n) can be any measurable Hermitian-valued function defined on the imaginary axis and the superscript ∗ denotes the complex conjugate transpose.Figure 2
The feedback system model structure.Lemma 2 (see [10]).
LetG(s)∈RH∞m×n and Δ be a bounded causal operator. If the following assumptions hold, (i)
for everyτ∈[0,1], the interconnection of G and τΔ is well-posed;(ii)
for everyτ∈[0,1], the IQC defined by Π is satisfied by τΔ;(iii)
there existsϵ>0 such that(11)GjωI∗ΠjωGjωI≤-ϵI,∀ω;then the feedback interconnection of G and Δ is stable.Remark 3.
From Lemma2, an important conclusion on IQC stability analysis can be drawn: if Δ satisfies several IQCs, it also satisfies all nonnegative linear combinations of IQCs. For example, if v, w satisfy the IQC defined by Π1 and the IQC defined by Π2, they will also satisfy the IQC defined by λ1Π1+λ2Π2 for all λ1,λ2≥0. In other words, the set of IQCs that Δ satisfies forms a “description” for the block Δ, and the more IQCs we know that Δ satisfies, the more precisely that we can describe the uncertainty of Δ.Lemma 4 (KYP lemma, [11]).
SupposeM(s)=Ccl(sI-Acl)-1Bcl+Dcl. Assume (Acl,Bcl) is stabilizable and Acl has no eigenvalues on the imaginary axis. Then the following conditions are equivalent: (i)
The systemM is stable, and M∞<γ.(ii)
We have∀ω∈[0,∞),(12)jω-Acl-1BclI∗CclTCclCclTDclDclTCclDclTDcl-γIjω-Acl-1BclI<0;(iii)
There exists a symmetric matrixP>0, such that(13)AclTP+PAclPBclBclTP-γI+CclTDclTCclDcl<0.Lemma 5 (see [11]).
By using linear fraction transformation, we convert the system configuration from Figure2 to Figure 3. Let e∈L2en[0,∞), Gl(s)∈RH∞ with corresponding dimension, and let the system L2-gain be the performance measurement.
Assume thatΔ satisfies an IQC defined by Π having the following block structure:(14)Π=Π11Π12Π12∗Π22;then the system in Figure 3 is stable and has robust L2-gain γ, if (i)
for everyτ∈[0,1], the interconnection of Gl and τΔ is well-posed;(ii)
for everyτ∈[0,1], the IQC defined by Π is satisfied by τΔ;(iii)
the frequency domain inequality,(15)GljωI∗ΓGljωI<0,holds for allω∈[0,∞), where (16)Γ=I0000Π11jω0Π12jω00-γ2I00Π12∗jω0Π22jω.Figure 3
The transferred system model structure.


## 3.2. Design Procedure

The algorithm for the design of Hill function considering both the stability and performance of the nonlinear system will be proposed.Proposition 6.
After shifting the equilibrium pointP (Figure 1) to the origin, Hill function h(x) takes the following form:(17)hx=βx+knkn+x+kn-β2,and there exist real scalars α1≥α2≥0, such that h(x) satisfies the IQCs defined by(18)π1=-2α1α2α1+α2α1+α2-2,π2=β22zy-y-z,where z≥0.Proof.
π
1 is from a sector bound condition α2v2≤h(v)v≤α1v2; to see this we notice that in the time domain(19)vhvTπ1vhv=2α1v-hvhv-α2v≥0.The derivation of π2 involves an approximation, in which h(v(t))=δ(v(t))v(t) with a scalar function δ∈L∞ with δ∞≤β/2; then(20)vhvTπ2vhv=vTβ22z-δ2z+δy-yv≥0.This ends the proof.Remark 7.
From Proposition6, we should know that the set of IQCs is a sufficient condition for stability; therefore, it is conservative. Yet the challenge for nonlinear system stability is significant, because all the existing methods, including Lyapunov theory and IQCs, are conservative.In Proposition6, π1 and π2 are for the case of scalar x and h(x). Here we consider the case of vectors x and h(x).Corollary 8.
Whenx=[x1,…,xi,…,xn,xn+1,…,x2n]T and H(x)=[0,…,0,h1(x1),…,hi(xi),…,hn(xn)]T, then IQCs for H(x) take the following forms:(21)Π1=Π1–11Π1–12Π1–21Π1–22,where (22)Π1–11=diag-2α11α21,…,-2α1iα2i,…,-2α1nα2n,Π1–12=diagα11+α21,…,α1i+α2i,…,α1n+α2n,Π1–21=Π1–12T,Π1–22=diag-2,…,-2,…,-2,Π2=Π2–11Π2–12Π2–21Π2–22,where (23)Π2–11=diagz1,…,zi,…,zn,Π2–12=diagy1,…,yi,…,yn,Π2–21=-Π2–12T,Π2–22=diag-z1,…,-zi,…,-zn,with zi≥0. Combining Π1 and Π2, we obtain the following IQCs for H(x):(24)Π=Λ1Π1–11Λ1Π1–12Λ1Π1–21Λ1Π1–22+Λ2Π2–11Λ2Π2–12Λ2Π2–21Λ2Π2–22,for any n×n nonnegative diagonal matrices Λ1 and Λ2.Example 9.
Taken=1, n1=3, k=1, and β=1 as an example; therefore,(25)hx=x+131+x+13-0.5.Assume that the evolution of x is within a limited a range [-1,1]; then a sector condition of α1=0.84, α2=0.35 with δ∈[-0.5,0.5] can form a good bound on h(x), which is shown in Figure 4.Figure 4
Hill functionh(x) and its bounds. In this picture, the solid line is the considered Hill function and dashed lines are the corresponding bounds.Proposition 10.
Consider the system configuration in Figure5, let e∈L2em[0,∞), Gl(s)∈RH∞ be the genetic network, and x(·)=[x1(·),…,xn(·),xn+1(·),…,x2n(·)]T and H(x)=[0,…,0,h1(x1),…,hi(xi),…,hn(xn)]T represent the designed feedback Hill function. We denote with ΔH the uncertainty derived from H(x), which satisfies the IQCs in (24).
Then the system in Figure5 is stable and has robust L2-gain γ, if (i)
for everyτ∈[0,1], the interconnection of Gl and τΔH is well-posed;(ii)
for everyτ∈[0,1], the IQC defined by Π is satisfied by τΔH;(iii)
the frequency domain inequality,(26)GljωI∗ΓGljωI<0,holds for allω∈[0,∞], where (27)Γ=I0000Π11jω0Π12jω00-γ2I00Π12∗jω0Π22jω,andΠijs take the form in (24).Figure 5
The system with Hill feedback model structure.Proof.
The result follows from applying the IQCs defined by (24) in Corollary 8.Proposition 11.
Let the Hill coefficientsHc be the set of proper values of (β,k) for a Hill function; given any (β,k)∈Hc and n, a Hill function can be uniquely defined and hence a set of IQCs in the form of (24) can be formulated. Then the design of Hill function feedback can be reformulated into the following linear matrix inequality (LMI) problem: (28)infΛ1,2>0,Hc,ωγ,such that (26) is satisfied.The frequency-dependent inequality in (26) can be resolved by using YALMIP (https://users.isy.liu.se/johanl/yalmip/). Since the linear matrix inequality (28) has to be solved simultaneously for all frequencies in the frequency domain, the basic strategy is to start with a small set of points and then check whether the solution satisfies all ω. If it does not, add one or more frequency points to the previous set and solve LMI again. Another more efficient way to solve the LMI is to transform (26) into the time domain by using the KYP lemma of Lemma 4 if the state-space realizations of Gl and Π have been obtained.


## 4. Biologically Inspired Example

In this section, we provide an example to illustrate the algorithm and design procedures developed in Section3. Consider the following GRN:(29)x˙t=Axt+But+Et,ut=Hxt,in which (30)A=-0.10000-0.500a10-0.100a20-0.5,B=0GB~0,G=0110,where a1, a2, and B~ are randomly generated in MATLAB and a1, a2, and the elements in matrix B~ follow the even distribution on the interval [0,1]. In this simulation, the generated data are a1=0.8, a2=0.1, (31)B~=0.70.90.11.The system gain of the original system G0=(sI-A)-1B is 94.2653. Although the system is originally stable, its performance is far from satisfactory. In order to further investigate the effects of noise e on u, the procedure of designing H(x) is slightly modified in this example. Specifically, it is an undesirable result if a small e causes a large u. Therefore, rather than minimizing H∞ norm from e to x, we minimize H∞ norm from e to [x,αu], where α is a constant. With Hc=(β,k)∣β∈[0.01,100],k∈[1,100] ([8]) and α=2, n=2, we solve the LMI problem stated in (28), resulting in the feedback Hill function with (β1,β2)=[10,100], (k1,k2)=[1,10]. As a result, H∞ norm of the closed-loop system is upper bounded by γ1=72.3593. Compared with γ0, a distinct improvement has been made for the system performance in terms of H∞ norm.


## 5. Conclusion

This paper tackles the problem of designing robust genetic circuit via constructing parameters of the Hill function. Targeting to minimize the noise effect in the genetic regulatory network, a design procedure is proposed. With the help of IQC approach, stability and performance of the designed circuit are guaranteed.



---

*Source: 5757943-2017-04-06.xml*