## Abstract

In terms of a piecewise affine system representation, which is a kind of hybrid system model, this article discusses a series of approaches to modelling, analysing and synthesizing a biological network such as a gene-regulatory network. First, the input assignment problem, the controllable state set problem (CSP) and the input trajectory generation problem are emphasized as control problems to be addressed for biological networks. Subsequently, after the modelling issue on biological networks developed in the systems and control community is briefly explained, the CSP is described in detail with reference to control of the quorum-sensing system in the pathogen *Pseudomonas aeruginosa*. Finally, an optimal control design method to the quorum-sensing system is proposed as a solution to the input trajectory generation problem.

## 1. Introduction

Since the mid-twentieth century, a piecewise affine (PWA) system-based approach to analysing system behaviours of biological networks has been highly appreciated (e.g. Glass & Kauffman 1973; Keener & Sneyd 2008). A typical case of this approach is that various kinds of findings are obtained by determining an analytical solution of the PWA system approximately representing them. However, the approach has been mostly applied to the class of biological systems expressed by autonomous differential equations with small dimension. On the other hand, in the last decade, the control theory of hybrid systems has been intensively developed in the systems and control community. In particular, for the class of PWA systems, one of the standard hybrid systems, various kinds of research issues on controllability, observability, stabilization, optimal control and model predictive control have been studied in addition to the analysis of autonomous systems such as stability.

Because of this situation, the research in applications of the above control theory to biological networks has brought in new developments of control theory suitable to biological networks. As for the identification/modelling issue in the systems and control community, PWA system identification methods for gene networks (Ferrari-Trecate *et al.* 2003; Drulhe *et al.* 2008) and PWA system approximation methods for nonlinear differential equations expressing a kind of chemical reaction such as Michaelis–Menten kinetics (Azuma *et al.* 2010) have been proposed. Reachability analysis and controllability analysis have also been developed by taking account of the aspects of biochemical reactions such as multi-affine functions and large uncertainties (Ghosh *et al.* 2003; Camlibel *et al.* 2007; Azuma *et al.* 2008; Batt *et al.* 2008). Furthermore, a kind of parameter-sensitivity analysis has been studied (Batt *et al.* 2007). In this way, the PWA systems approach is one of the most prospective methods in addressing control problems of biological networks.

In this article, we describe a series of PWA system model-based approaches to modelling, analysing and synthesizing biological networks. First, we discuss what kinds of problems should be addressed for biological networks from the systems and control points of view and introduce three kinds of fundamental problems, i.e. the input assignment problem (IAP), the controllable state set problem (CSP), and the input trajectory generation problem (ITP). This article mainly focuses on the second and third problems. The first problem has been treated by Sugiyama & Imura (2007). After explaining the modelling/identification issue via PWA system models in a brief way, we state a method for solving the CSP of biological networks, and its application to the quorum-sensing system in *Pseudomonas aeruginosa* as a case study. Finally, an optimal control design method to the quorum-sensing system is proposed as a solution to the ITP.

**Notation.** Let , and denote the real number field, the set of non-negative integers and the set of positive integers, respectively. We use the vector inequality *x*_{1}≤(<)*x*_{2} to express that each element of *x*_{1}−*x*_{2} is non-positive (negative). For the sets and , expresses the difference set. The set given as the form is called here the *polyhedron*, where *A*, *C* are some matrices and *b*, *d* are some column vectors.

## 2. Hybrid control theory for biological networks

First of all, we discuss what kinds of control problems should be addressed for analysing and controlling biological networks, such as gene-regulatory networks.

Consider the biological network expressed by the system defined on 2.1where is the state and is the control input, which are usually given by the concentration of chemical substances such as proteins. The notion of controllability of this system plays a fundamental role in many control problems such as stabilization by the state feedback and optimal control with the terminal state specified, usually defined as follows: for a given terminal state (or terminal set ), system (2.1) is said to be controllable at the initial state *x*_{0} at time *t*_{0} if for some *t*_{f}(>*t*_{0}) there exists an input trajectory *u*(*t*), *t*∈[*t*_{0},*t*_{f}] such that *x*(*t*_{f})=*x*_{f} (or ).

Thus, the controllability problem is defined as a kind of feasibility problem: given *x*_{0} and *x*_{f} (or ), determine if or not system (2.1) is controllable at *x*_{0}.

Indeed, the above problem is practical for a certain situation in considering the control of a gene network, but it will often be unsatisfactory for biologists and pharmacologists to use the solution as a clue or policy in estimating what experiment is significant.

First, for most biological networks, there are no input channels. Thus it will be useful to provide any information on which input channel is relatively effective. The problem of finding an effective input assignment should be addressed.

Next, most of the systems to be studied for biological networks are nonlinear. Furthermore, their mathematical models in general include large uncertainties on parameters/system structures and noise effects. These facts imply that even if we could determine that the model is controllable at a certain *x*_{0}, the real system may not. Thus, it will be of importance to determine a controllable state set, that is, a set of the initial state at which the model is controllable. From the shape and volume of the set, the ease of control as well as the feasibility of whether it is controllable at the current state will be able to be roughly estimated for the control channels in question.

Finally, it will be very useful to find an input trajectory for which the system is driven from the initial state *x*_{0} to the desired state *x*_{f} (or the desired state set ), and hopefully which is optimal in the sense that the least side-effect is produced. For example, qualitative findings on the optimal input trajectory may provide some clues for efficient experiments. In particular, this kind of information in the multi-input case will considerably decrease a huge number of experimental trials caused by the combinatorial explosion of input trajectories.

Thus, the following controllability problems will be fundamental for various control issues of biological networks.

— IAP: find an input assignment for the system to be controllable at every initial state in a set specified by a user.

— CSP: find a set of the controllable initial state.

— ITP: find an admissible input trajectory driven from a given initial state to a desired state or set.

The above problems in general involve state/input constraints. For example, a kind of spike (impulse) input of a substance as well as non-negativity of inputs may be imposed. Multiple injection of siRNA in a tumour model is an example of this case (Arciero *et al.* 2004). A dose of medicine is also considered as a discrete-valued input in some cases. The concentration of compounds, which is the state variable, is also non-negative. Furthermore, it may be required that some state variable should take a value in a prescribed bounded region so as to produce no side-effects on unmodelled dynamics.

Such state/input constraints as well as high nonlinear complexity usually render the above problems intractable. In addition, as mentioned before, the mathematical model of a biological network includes large uncertainties. Therefore, it may be a better strategy to approximately solve them rather than to rigorously solve them. A PWA system approach, which is discussed in the current article, will be one of the prospective tools in this line.

In §3, we briefly describe how we can derive PWA systems of biological networks.

## 3. Piecewise affine approximation of biological networks

Mainly two approaches to deriving PWA models of biological networks have been developed in the systems and control community: an identification approach and a function approximation approach. The former approach allows us to derive a PWA model directly from input–output data of the biological network as an identification problem (e.g. Ferrari-Trecate *et al.* 2003; Roll *et al.* 2004; Juloski *et al.* 2005; Nakada *et al.* 2005; Drulhe *et al.* 2008). Although many approaches have been developed, they in general include three steps: the input–output data-clustering based on *k*-means method or expectation maximization method, the estimation of sub-regions of the state space corresponding to each mode via support vector machine method, and the parameter identification of a linear (affine) system in each mode. However, it is still an open problem to systematically determine the number of modes as well as the dimension of the system.

On the other hand, in the latter approach, a PWA model is derived by approximating the nonlinear dynamics of biological networks often induced as chemical reactions such as Michaelis–Menten kinetics with PWA functions. This approach has been often used in analysing the solution behaviours of systems with small dimension. Recent developments on this issue provide various approaches (e.g. Henzinger *et al.* 1998; Girard 2002; Asarin *et al.* 2003). For example, a Lebesgue integral-based approach along this line has been proposed by Azuma *et al.* (2010), where the state space is partitioned according to the variation of the vector field as shown in figure 1. As a result, a PWA-approximated system with a smaller number of modes (i.e. partitions of the state space) can be obtained when compared with the case of evenly sized partition. Based on this notion, a design method of constructing a PWA approximation model with specified accuracy has been proposed.

By increasing the number of modes (i.e. affine systems), we can approximate the nonlinear dynamics with arbitrarily small discrepancy. However, the analysis of the resultant PWA model becomes much harder as the number of modes increases because the computational complexity exponentially increases with respect to the number of modes. Thus it is remarked that there is a trade-off to be made.

## 4. Controllability analysis of biological networks

This section reviews the CSP and its applications to the case of *P. aeruginosa*, which have been proposed by Azuma *et al.* (2008).

### (a) Controllable state set problem

Suppose that a discrete-time PWA system model (approximately) expressing the dynamics of a biological system is given:4.1where is the state, is the control input, is the discrete time, is the mode (), is the set of the mode values, is the number of the mode values and , , are constant matrices for mode *I*. In addition, denotes the sub-region of the state assigned to , given by the polyhedronfor given matrices and vectors , , and . For this sub-region, it is assumed that and for every , *I*≠*J* to guarantee that *I* is uniquely determined for each *x* and *u*. For the directed graph *Ξ* expressing for each current mode the allowable modes to which the mode is switched next, a set of the directed paths of length *T* from *I*_{0} and having at most switching times is denoted by .

Thus, the notion of a CSP is defined as follows.

### Definition 4.1

For system *Σ*, suppose that the controllability specification is given, where is the final time, is the initial state set to be examined, *Ξ*∈{0,1}^{M×M} is the directed graph expressing the mode switching rule, is the mode switch number, is the admissible state sequence set and is the admissible input sequence set. Then the set of all for which there exists an input vector sequencesatisfying and , where and , under the initial state *x*(0)=*x*_{0} is said to be an ‘S-controllable set’.

For example, if is defined as for the final state set , the S-controllable set represents the initial state set in which each state can reach the final state set at the *T*-th time with the satisfaction of the mode/input constraints. Thus, the CSP is formulated as the problem of finding an S-controllable set. However, it will be difficult in general to find it in a rigorous way owing to computational complexity. Thus, practically we may find a projection of the S-controllable set onto the state space of relatively small dimension, which we call here the ‘projected S-controllable set’, without calculating the whole S-controllable set.

Various kinds of problems in the controllability analysis of biological systems are covered by the above problem with . For example, the time-invariant state constraint , *k*∈{1,…,*T*} is expressed as . When an input is intermittently applied to a biological system as in the case of multiple injections of siRNA in the tumour model (Arciero *et al.* 2004), the input constraint at time *k*, where , is expressed by, for example, for every and otherwise , where is an input timing set given in advance. The volume of an S-controllable set is easily obtained as an index of controllability if , and are given by polyhedra; the corresponding S-controllable set is also given as a polyhedron.

Consider the parameter-driven controllability as a more specific example. One of the standard approaches to analysing behaviours of biological systems is based on numerical simulations. In this case, numerical simulations will be executed for a fixed initial state, but for various values of some physical parameters, to verify what parameter value can realize the prescribed solution behaviours of the system. On the other hand, if CSP is applied under *B*_{I}≡0 to such a situation, an initial state *set*, i.e. an S-controllable set, satisfying prescribed behaviours will be obtained for each parameter value; thus, the volume of the S-controllable set will be used as an index for determining what parameter value is effective.

Let us describe how to solve CSP very roughly, developed by Azuma & Imura (2007*a*,*b*). Since the problem is NP-hard, we use the randomized approach, where the mode trajectory is randomly selected. For each randomly selected , we first consider to calculate an S-controllable set represented aswhere matrices , , are given based on the state equation and the constraints. Noting thatwhere denotes the projection onto , andwe see that the problem is reduced into the projection problem of . This may be solved by well-known methods such as vertex enumeration and Fourier–Motzkin elimination for the problem with relatively small size of *n*, *m*, *T* or the linear programming-based method developed by Azuma & Imura (2007*a*,*b*), which provides a polynomial time algorithm with respect to *m* and *T*. Thus, we can calculate an approximate set of the -controllable set aswhere is a randomly selected subset of . The resultant set is interpreted in a certain probabilistic sense. It is remarked that the above calculation does not allow us to derive an explicit input trajectory, but characterizes the set of the initial state at which the input trajectory in question exists for a given mode trajectory.

### (b) Case study: quorum-sensing system in *Pseudomonas aeruginosa*

In this subsection, the application of the above method to the controllability analysis of the pathogen *P. aeruginosa* is explained.

#### (i) Quorum-sensing system model

The quorum-sensing system is a signal transduction mechanism in pathogens such as *P. aeruginosa*, *Serratia marcescens* and *Vibrio cholerae*, which works as a self-defence system for antibacterial factors in the environment such as plant and animal tissues. This kind of bacterium observes its own density through the concentration of self-produced, membrane-permeant signalling molecules (called the autoinducers) diffused in the extracellular space as well as in the intracellular space. Then if this density reaches a threshold level, the expression of several virulence genes is activated to produce bacterial exotoxins and/or to form biofilms.

This article introduces an application to *P. aeruginosa*, which is a major opportunistic human pathogen living in various environments including human bodies, seawater, soil, and so on. It is one of the crucial issues to find an effective treatment in spite of the resistance to antibiotics and disinfectants caused by the quorum-sensing system in this pathogen. The gene-regulatory network of the quorum-sensing system in *P. aeruginosa* is shown in figure 2 (Fagerlind *et al.* 2003, 2005). This network consists of two subsystems: the *las* system and the *rhl* system. In the *las* system, the proteins LasR, LasI and RsaL are coded by the genes *lasR* (transcriptional activator), *lasI* (autoinducer synthase) and *rsaL* (LasI inhibitor), respectively. If the concentration of the autoinducer OdDHL reaches a threshold level, the LasR/OdDHL complex is produced to induce the expression of several virulence genes as well as *lasI* and *rsaL*, while *rsaL* inhibits the expression of *LasI*. The *rhl* system composed of the genes *rhlR* and *rhlI* also works in a similar way to the *las* system. The LasR/OdDHL complex in the *las* system induces the expression of *rhlR*, and OdDHL works as an *rhl* system inhibitor.

The following mathematical model of this system has been recently developed by Fagerlind *et al.* (2003):4.2where *x*_{c1}, *x*_{c2}, *x*_{c3}, *x*_{r1}, *x*_{r2}, *x*_{s}, *x*_{a1} and *x*_{a2} are the cellular concentrations of the LasR/OdDHL complex, the RhlR/BHL complex, the RhlR/OdDHL complex, LasR, RhlR, RsaL, OdDHL and BHL, respectively, whose units are [μM]. The model has three equilibrium states for *x*:=[*x*_{c1}*x*_{c2}*x*_{c3}*x*_{r1}*x*_{r2}*x*_{s}*x*_{a1}*x*_{a2}]^{⊤}, where two states are asymptotically stable and the other is unstable. One of the two stable equilibrium states corresponds to the case that the virulence gene is activated, while the other, denoted by *x*_{eu}, is the case that is not activated.

#### (ii) Controllability problem

The use of an OdDHL antagonist such as furanone C30, which inhibits the binding of LasR and OdDHL, has been currently considered as one of the most effective methods to prevent the expression of virulence genes via the quorum-sensing system (Givskov *et al.* 1996; Wu *et al.* 2004). A specific antibody for capturing some of the extracellular OdDHL is also used to reduce the intracellular concentration of OdDHL via a diffusion effect.

Which gene expressions will be the most essential for controlling the expression of virulence genes? To give an answer, by means of solving CSP, we consider here the problem of finding an input channel at which there exists a control input function such that the virulence genes are not activated, i.e. the state of *P. aeruginosa* is driven to a set around the second equilibrium point *x*_{eu} in a finite time.

It turns out that the nonlinear model (4.2) can be expressed by the 8-modal discrete-time PWA system model of (4.1) approximately with a high degree of accuracy, where *I*∈{1,2,…,8}, the sub-regions and the system matrices , are determined based on the similarity of the vector fields and of the solution behaviours by numerical simulations (see Azuma *et al.* (2008) for more details). We have four cases of the input term *B*^{c}*u*, which corresponds to applying some antagonist or specific antibody acting on the deviation of the state specified by *B*^{c}, i.e.

— case 1:

*B*^{c}:=[1 0 0 0 0 0 0 0]^{⊤}(for LasR/OdDHL complex),— case 2:

*B*^{c}:=[0 0 0 1 0 0 0 0]^{⊤}(for LasR),— case 3:

*B*^{c}:=[0 0 0 0 0 0 1 0]^{⊤}(for OdDHL), and— case 4: (for both LasR and on OdDHL).

We also suppose , *k*=0,1,…,*T*−1, for each input (note that the antagonists/specific antibodies are a kind of depressor). Then for each of the above cases, we are interested in the S-controllable sets with the following : *T*:=20, ,*σ*=8, for and () and for defined above. Here, *Ξ* and *σ* are given so as to remove obviously unnatural behaviours. Under the above conditions, the problem is solved for the (*x*_{c1},*x*_{r1},*x*_{a1})-space.

The resultant projected S-controllable sets are shown in figure 3. Because of the use of the randomized algorithm, the set except for the polyhedra there expresses the projection of the set of satisfying , where . All results have been obtained within 1–4 h, where MATLAB 7.1.0 and MPT toolbox 2.6.1 (including CPLEX solver for linear programming) are used with 64-bit Xeon (dual core) × 4 (3.66 GHz) and 8 GB memory (DDR2/400). From the above results, we see that case 1 is the best among the cases examined here. Biologists also consider that it is important to reduce the concentration of this complex (Fagerlind *et al.* 2003, 2005). Furthermore, from the results of cases 2 and 3, we see that it is difficult to inhibit the expression of virulence genes by controlling the concentration of only one of LasR and OdDHL. Thus, although the OdDHL-antagonist/specific antibody is a major approach in the field of biology, this result suggests that the claim that the OdDHL-antagonist/specific antibody is effective (Fagerlind *et al.* 2003, 2005) may not be always true. More precisely, we conclude that, under the input constraint , the input of cases 2 and 3 will be effective when the cellular concentration of the LasR/OdDHL complex is small. Finally, it turns out from case 4 that the simultaneous use of two inputs is effective similar to that of case 1. The above results show that an antagonist/specific antibody for one of two proteins from which the complex is produced may not be enough to reduce the concentration of the complex.

## 5. Optimal control of quorum-sensing system

In this section, we illustrate an example of optimal control input design for the *las* system of *P. aeruginosa* depicted in figure 2. We ignore here the *rhl* system for brevity.

### (a) Piecewise affine model

We begin with some experimental data obtained in the following two different situations: (i) with *rsaL* and (ii) without *rsaL* (this is done by making *rsaL* defective). Figures 4 and 5 show the concentration of the autoinducer and the scaled population of *P. aeruginosa*.

We move to mathematical modelling of this quorum-sensing mechanism. Let us introduce the variables *x*_{0} (scaled population of *P. aeruginosa*), *x*_{a} (concentration of the autoinducer OdDHL), *x*_{c} (cellular concentration of the LasR/OdDHL complex), *x*_{r} (concentration of LasR) and *x*_{s} (concentration of RsaL). Then the following nonlinear mathematical model can be obtained:
5.1
5.2
5.3
5.4
5.5
where *u* is the control input that inhibits the expression of autoinducer, and *M*_{x0} denotes the stationary value of *x*_{0}. Figures 6 and 7 show the integral curves of *x*_{a} and *x*_{0} of this model. As in the experimental data, we performed two cases: (i) with *rsaL* (*x*_{s}(*t*)>0) and (ii) without *rsaL* (*x*_{s}(*t*)≡0). We can see that this model emulates the behaviour of the experimental data in both cases. It is remarked that the above model is slightly different from the *las* system model of (4.2) in the sense that the population of *P. aeruginosa* is time-varying in the above model, not stationary.

In a similar way to discussions in the previous sections, we perform PWA approximations. We first replace equations (5.2)–(5.5) by a continuous-time PWA system model. For simplicity, the state space is partitioned into 23 regions according to the value of only *x*_{a} and *x*_{c}, i.e.with and appropriately defined. System matrices are determined by the Taylor expansion of the nonlinearity around a certain point in . Then we discretized the PWA system with respect to time with the sampling time 1 (h) to derive a discrete-time PWA system model of (4.1). Since *x*_{0}(*t*) in equation (5.1) evolves autonomously, we approximate *x*_{0}(*t*) by a step function as in figure 8. Figure 9 shows the integral curves of *x*_{a} and *x*_{c} of the original nonlinear system of equations (5.2)–(5.5) and the approximated discrete-time PWA system for a certain initial state and *u*≡0.

### (b) Design example

Our goal is to inhibit the increase of *x*_{c} by suitably adding the control input *u*. In view of this, we search positive step signals of period 1 (h) for *k*=0–10 for control input that drives the initial state [*x*_{0} *x*_{a}*x*_{c}*x*_{r}*x*_{s}]=[0.05 0.05 0.05 0.05 0.05] in [0,0.09]^{5} at *k*=10. To minimize the amount of the input, we employ the cost function for the optimization. To solve this problem, we use the mixed logical dynamical (MLD) model (Bemporad & Morari 1999):andwhere is the state, , is the control input, and are auxiliary continuous and binary variables, respectively, and system matrices are time-varying. The discrete-time PWA system model of (4.1) can be equivalently expressed by the above MLD model. See Bemporad & Morari (1999) for further details. Note that the above MLD model has linear equations and linear constraints, although it includes binary variables. Thus the problem of minimizing the cost function *J* under the PWA system with respect to the control input is reduced into the mixed integer linear programming (MILP) problem, which can be solved by an efficient solver such as ILOG CPLEX solver.

The obtained input *u*(*k*) and the resultant *x*_{c}(*k*) are shown in figures 10 and 11, respectively. We can see that *x*_{c} is suppressed by adding input as desired (cf. *x*_{c} in figure 9). The shape of the input sequence is slightly counterintuitive in that it has an unnatural peak at *k*=6. We attempt to give some interpretation for this. In equation (5.2), the input affects only through *x*_{a}*u*. In view of figures 7 and 12, we see that the control input sequence is designed so as to smoothly increase the effect of *x*_{a}*u* as the population of *P. aeruginosa* increases around *k*=6. In this way, it turns out that the least amount of control input should be obtained heavily depending on the state *x*_{0} and *x*_{a}.

## 6. Conclusion

This article reviewed a PWA system model approach to analyse and synthesize biological networks such as gene-regulatory networks from the viewpoints of the systems and control community. The emphases are on controllability and optimal control issues for biological networks, whose solutions can be given by exploiting a PWA system model instead of an original nonlinear model, and mixed integer mathematical programming. This kind of approach is still ongoing and there are many significant unresolved issues to be addressed towards the developments of control theory for biological networks including various kinds of applications.

## Acknowledgements

The first author would like to express his thanks to Prof. Kazuyuki Aihara, University of Tokyo, for his fruitful comments. This research is partially supported by the Japan Society for the Promotion of Science (JSPS) through its Funding Programme for World-Leading Innovative R&D on Science and Technology (FIRST Programme).

## Footnotes

One contribution of 12 to a Theme Issue ‘Theory of hybrid dynamical systems and its applications to biological and medical systems’.

- © 2010 The Royal Society