# Time-optimal control of piecewise affine bistable gene-regulatory networks: preliminary results

Nicolas Augier, Agustín Gabriel Yabo

► To cite this version:

Nicolas Augier, Agustín Gabriel Yabo. Time-optimal control of piecewise affine bistable gene-regulatory networks: preliminary results. ADHS 2021 - 7th IFAC Conference on Analysis and Design of Hybrid Systems, IFAC, Jul 2021, Brussels, Belgium. 10.1016/j.ifacol.2021.08.499. hal-03099681v4

HAL Id: hal-03099681

https://inria.hal.science/hal-03099681v4

Submitted on 9 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Time-optimal control of piecewise affine bistable gene-regulatory networks: preliminary results*

Nicolas Augier* Agustín Gabriel Yabo*

* Université Côte d'Azur, Inria, INRAE, CNRS, Sorbonne Université, Biocore Team, Sophia Antipolis, France. nicolas.augier@inria.fr, agustín.yabo@inria.fr.

**Abstract:** In this preliminary work we give a geometric characterization of the time-optimal trajectories for a piecewise affine bistable switch. Such hybrid models play a major role in systems biology, as they can expressively account for the behaviors of simple gene-regulatory networks. A key example of these biological circuits is the genetic toggle switch, for which the optimality of a state switch has not been studied to the present day. The main result of this work to prove that the optimal control strategy is a simple bang-bang feedback control law. In opposition to previous works, we show that a time-optimal transition between states should pass by an undifferentiated state, well known in cell biology for playing a major role in fate differentiation of cells. We conclude by showing numerical simulations of optimal trajectories illustrating the structure of the bang-bang optimal control for different scenarios.

**Keywords:** Genetic Regulatory Systems; Hybrid systems; Biological systems; Optimal control; Genetic toggle switch

## 1. INTRODUCTION

Controlling systems of monotone feedback loops (defined for instance in Glass and Pasternack (1978); Mallet-Paret and Smith (1990)) is a critical task for many applications in biology (Chaves and Gouzé, 2011; Chambon et al., 2020), notably in Synthetic Biology, in the understanding of gene-regulatory mechanisms, and for pharmacological treatments. Such systems involve several variables $(x_i)_{i \in \{1, \dots, n\}}$ representing gene expressions, which interact through the general dynamics $\dot{x}_i = f_i(x_i) + g_i(x_{i-1})$, $i \in \{1, \dots, n\}$, with periodic conditions $x_0 = x_n$, where $(g_i)_{i \in \{1, \dots, n\}}$ are monotone functions from $\mathbb{R}$ to $\mathbb{R}$. The interaction between $x_{i-1}$ and $x_i$ is said to be positive (respectively, negative) when $g_i$ is increasing (respectively, decreasing). These systems may present "typical" behaviours, which have been studied in Mallet-Paret and Smith (1990); Farcot and Gouzé (2008); Chaves and Gouzé (2011); Glass and Pasternack (1978); Farcot and Gouzé (2008); Chambon et al. (2020), which are essentially the multistability or the presence of a limit cycle. Controlling such systems by a real non-negative control input $u$ acting on the interaction terms has been initiated in particular frameworks in Chaves and Gouzé (2011); Chambon et al. (2020), and in both cases the controlled dynamics writes $\dot{x}_i = f_i(x_i) + g_i(x_{i-1}, u)$, $i \in \{1, \dots, n\}$, with periodic conditions $x_0 = x_n$ and $(g_i)_{i \in \{1, \dots, n\}}$ being functions from $\mathbb{R}^2$ to $\mathbb{R}$, which are monotone w.r.t. the first variable, for every $u \ge 0$. Notably, in Chambon et al. (2020),

the goal was to suppress undesirable oscillations, and in Chaves and Gouzé (2011) a control strategy was proposed which allowed to achieve transitions between the steady states of a two dimensional bistable piecewise affine system. Here, we focus on the model from Chaves and Gouzé (2011), which consists in a piecewise affine classical model describing a bistable switch. In synthetic biology, such systems are widely used to represent synthetic gene-regulatory networks which, in practice, are capable of forming biological memory units. These synthetic circuits, known in the literature as genetic toggle switches, have been first implemented experimentally by Gardner et al. (2000) in E. coli through two genes lacI and tetR mutually repressing each other. In this context, understanding how to regulate them has become highly relevant for their vast implications in biotechnology and biocomputing. In the mathematical modeling framework, functions $(g_i)_{i \in \{1,2\}}$ are considered as step functions with fixed thresholds corresponding to negative interactions, and $f_i(x) = -\gamma_i x$, where $\gamma_i > 0$ for every $i \in \{1, \dots, n\}$. $g_1(x) = k_1 s^-(x, \theta_2)$, $g_2(x) = k_2 s^-(x, \theta_1)$, $k_j > 0$ for every $j \in \{1, 2\}$, where for $\theta \in \mathbb{R}$, $s^-(\cdot, \theta) : \mathbb{R} \to \mathbb{R}$ is such that

$$s^-(x, \theta) = \begin{cases} 1 & \text{if } x < \theta, \\ 0 & \text{if } x > \theta. \end{cases}$$

The authors of Chaves and Gouzé (2011) insisted on the role of a particular state, which can be considered, biologically, as an "undifferentiated state", where no gene is predominant, and from which the system can evolve towards one of the two attractors. From a biological point of view, such state seems to play a key role in cell decision making and cell fate differentiation (Balázsi et al., 2011). In recent works, Lugagne et al. (2017) constructed a real-time control scheme that allowed to stabilize a genetic

<footnote>* This work was partially supported by ANR project Maximic (ANR-17-CE40-0024-01), Inria IPL Cosy and Labex SIGNALIFE (ANR-11-LABX-0028-01). We also acknowledge the support of the FMJH Program PGMO and the support to this program from EDF-THALES-ORANGE.</footnote>
trigger switch around this undifferentiated point. This was done by externally adding the diffusible molecules IPTG<fnref n="1" /> and aTc<fnref n="2" />, which are known to repress the lacI and tetR genes. Mathematically, the instability of the undifferentiated state is modeled in Chaves and Gouzé (2011) through the presence of a Filippov non-smooth "saddle" singularity (Filippov, 1960).

In this paper, we investigate the time-optimal control strategies for the aforementioned bistable system. Our aim is to induce transitions between the two steady states in minimal time. While most of the work in the subject has been dedicated to the feasibility and implementation of such systems, the optimality aspect has received little or no attention from the community. Time-optimal transfers between two steady states of a bistable system has been studied in other settings, for instance for two-level quantum systems in Boscain and Mason (2006): Boscain et al. (2014), where the authors aimed at inducing efficient transitions between two quantum states. In our case, the steady states of the system cannot be reached in finite time, hence one has to consider "partial targets", that is, driving a given protein to a certain fixed value larger than its corresponding threshold. In this regard, we show that time-optimal strategies for such a problem have a very specific geometric description. When the initial state is far enough from the target, that is below a curve called separatrix, they consist in a concatenation of two bang arcs, and the optimal trajectories follow:

* a first phase in which the system reaches the separatrix;

* a second phase where the system slides along this curve, until reaching the "undifferentiated" point of the biological system in finite time;

* a third phase, where the system leaves this curve, slides along a second fixed curve and reaches its target.

These two curves correspond to the stable and unstable manifolds of the undifferentiated "saddle type" singularity, and the point where the dynamics achieves its transfer is nothing but the corresponding Filippov equilibrium. The latter behavior can be compared to the turnpike phenomenon (Trélat and Zuazua, 2015), where the optimal trajectory for a given OCP (Optimal Control Problem) for large final times is shown to remain close to a steady-state trajectory solution of the associated static OCP.

## 2. BISTABLE-SWITCH MODEL

### 2.1 Free dynamics

Consider two variables $x_1$ and $x_2$ which represent two genes mutually inhibiting each other. The individual dynamics, defined in Filippov sense (see Chaves and Gouzé (2011, Appendix A) for more precisions), is the following

$$ \begin{cases} \dot{x}_1 = -\gamma_1 x_1 + k_1 s^-(x_2, \theta_2), \\ \dot{x}_2 = -\gamma_2 x_2 + k_2 s^-(x_1, \theta_1), \end{cases} \quad (1) $$

where for $j \in \{1, 2\}$, $x_j \in \mathbb{R}$, and for $\theta \in \mathbb{R}$, $s^-(\cdot, \theta) : \mathbb{R} \to \mathbb{R}$ is such that

$$ s^-(x, \theta) = \begin{cases} 1 & \text{if } x < \theta, \\ 0 & \text{if } x > \theta. \end{cases} $$

It is assumed that $s^-(x) \in [0, 1]$ for $x = \theta$. The positive constants $(\gamma_j)_{j \in \{1, 2\}}$, $(k_j)_{j \in \{1, 2\}}$ correspond, respectively, to the degradation and the production rates of each variable. It is a classical fact (see Chaves and Gouzé (2011)) that the domain $K = [0, \frac{k_1}{\gamma_1}] \times [0, \frac{k_2}{\gamma_2}]$ is forward invariant by the dynamics of Equation (1). From now on we consider only solutions evolving in $K$. Define the regular domains

$$ \begin{align*} B_{00} &= \{(x_1, x_2) \in \mathbb{R}^2 \mid 0 < x_1 < \theta_1, 0 < x_2 < \theta_2\}, \\ B_{01} &= \{(x_1, x_2) \in \mathbb{R}^2 \mid 0 < x_1 < \theta_1, \theta_2 < x_2 < \frac{k_2}{\gamma_2}\}, \\ B_{10} &= \{(x_1, x_2) \in \mathbb{R}^2 \mid \theta_1 < x_1 < \frac{k_1}{\gamma_1}, 0 < x_2 < \theta_2\}, \\ B_{11} &= \{(x_1, x_2) \in \mathbb{R}^2 \mid \theta_1 < x_1 < \frac{k_1}{\gamma_1}, \theta_2 < x_2 < \frac{k_2}{\gamma_2}\}. \end{align*} $$

Equation (1) restricted to a regular domain is an affine dynamical system on $\mathbb{R}^2$ having an asymptotically stable equilibrium, called focal point for System (1). Each region $B_{ij}$ for $i, j \in \{0, 1\}$ has a focal point $\phi_{ij} = (\bar{x}_i, \bar{x}_j)$ corresponding to

$$ \bar{x}_i = \frac{k_i}{\gamma_i} s^-(\bar{x}_j, \theta_j). $$

Thus, system (1) has two locally asymptotically stable steady states

$$ \phi_{10} = \left( \frac{k_1}{\gamma_1}, 0 \right) \in B_{10}, \quad \phi_{01} = \left( 0, \frac{k_2}{\gamma_2} \right) \in B_{01}, $$

and an unstable Filippov equilibrium point at $(\theta_1, \theta_2)$.

<figure><img src="image_1.png" /><figcaption>Fig. 1. Stream plot with free dynamics. System parameters are $\gamma_1 = 1.2$, $\gamma_2 = 1.8$, $\theta_1 = 0.6$, $\theta_2 = 0.4$, $k_1 = k_2 = 1$.</figcaption></figure>

The hybrid nature of system (1) is a result of the step function $s^-$, which allows a qualitative analysis of the classical non-hybrid bistable switch of Gardner et al. (2000).

### 2.2 Controlled dynamics

In accordance with experimental studies (Gardner et al., 2000), we write the controlled dynamics assuming that the

<footnote n="1">¹ isopropyl-β-D-thiogalactopyranoside</footnote>

<footnote n="2">² anhydrotetracycline</footnote>
synthesis rates of each gene can be externally catalyzed or inhibited. Indeed, it is rather common in experimental biology to affect the synthesis rate of messenger RNA by introduction of a plasmid (which are essentially small circular DNA molecules that can be constructed to include an inducible promoter of the studied gene). This way, the transcription rate can be externally modified by aggregation of an inducer, which we model through an external control multiplying the production rate of the variables. We consider the controlled system, defined in Filippov sense,

$$ \begin{cases} \dot{x}_1 = -\gamma_1 x_1 + u(t) k_1 s^{-}(x_2, \theta_2), \\ \dot{x}_2 = -\gamma_2 x_2 + u(t) k_2 s^{-}(x_1, \theta_1), \end{cases} \quad (S) $$

where the control $u(\cdot) \in L^\infty([0, t_f], [u_{\min}, u_{\max}])$, with $0 < u_{\min} < 1 \le u_{\max}$. We make the following assumptions on the parameters of the system (for more details, see Chaves and Gouzé (2011)).

**Assumption 1.** The parameters $(\gamma_j)_j$ and $(k_j)_j$ satisfy

$$ \theta_j < \frac{k_j}{\gamma_j}, \quad j \in \{1, 2\}; \quad \frac{\theta_2}{\theta_1} > \frac{k_2}{k_1} \frac{\gamma_1}{\gamma_2}; \quad \frac{\theta_2}{\theta_1} < \frac{k_2}{k_1}. $$

Note that Assumption 1 implies $\gamma_1 \neq \gamma_2$. We define the separatrix of a system, which plays a fundamental role in the global dynamics of both the open loop system (1) and controlled system (S).

**Separatrix** For a fixed value of $u(t) \equiv u \in [u_{\min}, u_{\max}]$, the separatrix ($S_u$), is defined as the stable manifold of the Filippov equilibrium ($\theta_1, \theta_2$) for Equation (S) restricted to the closure $\overline{B_{00}}$ of $B_{00}$. In the coordinates $(x_1, x_2) \in B_{00}$, for $u \ge 0$, the separatrix ($S_u$) can be written as the curve of equation

$$ \begin{aligned} x_2 &= \alpha(x_1, u) = \frac{k_2 u}{\gamma_2} - \left( \frac{k_2 u}{\gamma_2} - \theta_2 \right) \left( \frac{\frac{k_1 u}{\gamma_1} - x_1}{\frac{k_1 u}{\gamma_1} - \theta_1} \right)^{\frac{\gamma_2}{\gamma_1}} \\ x_2 &= \frac{k_2}{\gamma_2} - \left( \frac{k_2}{\gamma_2} - \theta_2 \right) \left( \frac{\frac{k_1}{\gamma_1} - x_1}{\frac{k_1}{\gamma_1} - \theta_1} \right)^{\frac{\gamma_2}{\gamma_1}}. \end{aligned} $$

($S_1$) Using the latter, we define the regions

$$ (S_u)^+ = \{(x_1, x_2) \in \mathbb{R}^2 \mid 0 < x_1 < \theta_1, \alpha(x_1, u) < x_2\}, $$

$$ (S_u)^- = \{(x_1, x_2) \in \mathbb{R}^2 \mid 0 < x_2 < \theta_2, \alpha(x_1, u) > x_2\}. $$

such that $K$ is divided into

$$ K = (\overline{S_u})^+ \cup (\overline{S_u})^- \cup \overline{B_{11}}. $$

The solutions of Equation (S) having initial conditions in $(S_u)^-$ (respectively, $(S_u)^+$) reach $B_{10}$ (respectively, $B_{01}$) in finite time. Moreover, $B_{10}$ (respectively, $B_{01}$) is included in the bassin of attraction of $\phi_{10}$ (respectively, $\phi_{01}$). Notice that, for a fixed value of $u(t) \equiv u \in [u_{\min}, u_{\max}]$, the solutions of Equation (S) having initial conditions in $(S_u)$ reach $(\theta_1, \theta_2)$ in finite time. Once having reached such a Filippov point, the solution of Equation (S) is then defined by differential inclusion in the Filippov sense. Roughly speaking, there exist several solutions that will reach either $B_{01}$ or $B_{10}$ (see (Chaves and Gouzé, 2011, Appendix) for more precise informations about Filippov solutions of such system).

# 3. TIME-OPTIMAL TRANSFER

## 3.1 Problem formulation

The genetic toggle switch represents the lowest level of granularity in more complex genetic circuits. The state of this device is determined by gene expression in the boolean form (low, high) and (high, low). The problem investigated in this work is to achieve a transition from one boolean state to the other in minimal time. In the mathematical context, the latter translates into finding trajectories that drive the solution $(x_1(t), x_2(t))$ of Equation (S) towards the steady states $\phi_{01}$ and $\phi_{10}$ of Equation (1) in minimum time (where these states correspond to the differentiated states aforementioned). However, due to the lack of controllability in direction $x_1$ (respectively, $x_2$) of Equation (S) restricted to $B_{01}$ (respectively, $B_{10}$), one has to relax the problem. More precisely, the steady state $\phi_{01}$ (respectively, $\phi_{10}$) cannot be reached in finite time, because $u$ does not act on $x_1$ in the domain $B_{01}$ (respectively, $x_2$ in the domain $B_{10}$). Thus, we will be first interested in driving $x_2(t)$ towards an arbitrary value $x_2(t_f) = x_2^f > \theta_2$ (for instance, the value $x_2^f = \frac{k_2}{\gamma_2}$ corresponding to the $x_2$-component of the steady state $\phi_{01}$), with the constraint that at the final time, $x_1(t_f)$ belongs to the interval $[0, \theta_1)$. This target choice ensures that, at the final time, the gene $x_2$ is strongly expressed while the gene $x_1$ is weakly expressed. The symmetric problem, which is equivalent, consists in driving $x_1(t)$ towards an arbitrary value $x_1(t_f) = x_1^f > \theta_1$, with the constraint that at the final time, $x_2(t_f)$ belongs to the interval $[0, \theta_2)$. Fix $x_1^0 \ge \theta_1$, $x_2^0 \le \theta_2$, $x_2^f > \theta_2$, and consider the minimization problem

$$ \left\{ \begin{array}{l} \text{minimize } t_f \ge 0 \\ x(t) = (x_1(t), x_2(t)) \text{ is subject to (S),} \\ x(0) = (x_1^0, x_2^0), \\ x_2(t_f) = x_2^f, \\ x_1(t_f) \in [0, \theta_1), \\ u(\cdot) \in [u_{\min}, u_{\max}]. \end{array} \right. \qquad (OCP) $$

## 3.2 Reachability of the terminal state

A fundamental aspect of OCPs with fixed terminal state is the existence of a solution. Such matter is directly linked to the reachability and controllability analysis of the dynamical system, which are often hard to conduct analytically. In this work, we provide sufficient conditions for the feasibility of the proposed trajectory, and we show that the piecewise constant control strategy presented in Chaves and Gouzé (2011) achieves the objective, serving as a candidate to (OCP). Such strategy drives the system from an initial state in $B_{01}$ to $\phi_{10}$ (or $B_{10}$ to $\phi_{01}$ for the symmetric problem). One can show that, under Assumption 1, there exist $u_{\min} < \frac{\theta_1 \gamma_1}{k_1}$ and $u_{\max} \ge 1$ such that

$$ \Phi^* \in (S_{u_{\max}})^+, $$

with

$$ \Phi^* \doteq \left( \frac{u_{\min} k_1}{\gamma_1}, \frac{u_{\min} k_2}{\gamma_2} \right). $$
Concerning the symmetric problem, one can show the existence of another choice of $u_{min} < \frac{\theta_1\gamma_1}{k_1}$ and $u_{max} \ge 1$ such that $\Phi^* \in (S_{u_{max}})^-$. For this paper, let us focus on the first case, and state the following assumption.

**Assumption 2.** Bounds $u_{min}$ and $u_{max}$ are chosen such that $\Phi^* \in (S_{u_{max}})^+$.

Based on the solution developed in Chaves and Gouzé (2011), we first propose an input control constrained to two possible values $\{u_{min}, u_{max}\}$ corresponding to the low and high synthesis control. The control law is expressed in terms of the state and time as

$$ u(x,t) = \begin{cases} u_{min} & \text{if } x \in B_{10}, \\ u_{min} & \text{if } t \in [0, t_s), x \in B_{00}, \\ u_{max} & \text{if } t \in [t_s, \infty), x \in B_{00}, \\ u_{max} & \text{if } x \in B_{01}. \end{cases} \quad (2) $$

for $T_1$ sufficiently large. During the first phase with $u \equiv u_{min}$, every focal point of the system belongs to $B_{00}$, hence the solution $x(t)$ of Equation (S) converges towards the point $\Phi^* \in B_{00}$ when $t \to \infty$. During the second phase with $u \equiv u_{max}$, state $x(t)$ reaches $B_{01}$ in finite time, and $x_2(t)$ converges towards $x_2^f$ in finite time. From that point, an open-loop control $u \equiv 1$ would drive $x(t)$ to $\phi_{01}$ when $t \to \infty$. Indeed, under Assumptions 1 and 2, and by choosing $T_1$ sufficiently large, the control strategy (2) ensures that any trajectory starting from $(x_1^0, x_2^0)$ reaches a final point meeting $x_1 \in [0, \theta_1)$ and $x_2 = x_2^f$ in finite time, which shows that the set of admissible controllers for problem (OCP) is non empty.

### 3.3 Hybrid optimal control problem with a fixed domain sequence

Consider two compact subsets $M_0$ and $M_1$ of $\mathbb{R}^2$, and assume $M_1$ is reachable from $M_0$ for system (S), that is, such that there exists a time $t_f > 0$, a control $u(\cdot) \in L^\infty([0, t_f], \Omega)$ and $x_0 \in M_0$ such that the solution $x(t)$ of Equation (S), defined in the Filippov sense with initial condition $x(0) = x_0$ satisfies $x(t_f) \in M_1$. Consider the problem of steering the system (S) from $M_0$ to $M_1$ in minimal time $t_f$. The latter can be represented by the following OCP subject to piecewise linear dynamics

$$ \left\{ \begin{array}{l} \text{minimize } t_f \\ x(t) \text{ is subject to (S),} \\ x(0) \in M_0, \\ x(t_f) \in M_1, \\ u(\cdot) \in \Omega. \end{array} \right. \qquad (\text{PWLOCP}) $$

In order to properly define the problem, one has to choose a sequence $\mathbf{P}$ in the set $\{B_{00}, B_{01}, B_{10}, B_{11}\}$ of regular domains.

**Definition 1.** Let $B = (B_j)_{j \in \{1, \dots, k\}}$ be a sequence of regular domains. We say that a solution $x(t)$ of Equation (S) is B-admissible if there exists a time $t_f > 0$, a control $u(\cdot) \in L^\infty([0, t_f], \Omega)$, and times $t_0 = 0 < t_1 < \dots < t_k \le t_f$ such that $x(t) \in B_j$ for every $t \in \Delta_j$, where $\Delta_j = (t_j, t_{j+1})$.

In particular, the previous definition excludes sliding modes along the frontier between two successive regular domains. Additionally, we require two more assumptions

related to the reachability of B-admissible solutions for the general case.

**Assumption 3.** $M_0$ (respectively, $M_1$) is included in the adherence of a regular domain $B_{jk}$ (respectively, $B_{qi}$), for $j, k, q, i \in \{0, 1\}$.

**Assumption 4.** Assume that there exists a time $T > 0$ and a B-admissible solution $(x(t), u(t))$ such that $x(0) \in M_0$ and $x(T) \in M_1$.

Notice that for given sets $M_0, M_1$, the choice of the sequence $B$ is not unique in general. Assume that $M_0, M_1, B$ satisfy the assumptions 3 and 4. For a fixed sequence $B = (B_j)_{j \in \{1, \dots, k\}}$, we can consider Problem (PWLOCP) restricted to B-admissible trajectories, and necessary conditions of optimality for this problem can be directly derived from the HMP (Hybrid Maximum Principle) stated in Theorem 3 of Appendix A.

## 4. MAIN RESULTS

We are interested in solving (OCP) among continuous B-admissible trajectories, as defined in Section 3.3, with $B = (B_{10}, B_{00}, B_{01})$ (which has been shown to be an admissible sequence of domains in Section 3.2), $M_0$ restricted to a point in $K$, and $M_1 = \{(x_1, x_2) \in K \mid x_1 \in [0, \theta_1), x_2 = x_2^f\}$ which has already been proven to be reachable in finite time, and so the assumptions 3 and 4 are satisfied. As previously said, the problem can be further analyzed by applying HMP. The Maximum Principle in the Hybrid framework requires to define functions $(\phi_i)_i$ and $(\eta_j)_j$ that guarantee the continuity of the trajectories and the changes of dynamics at the frontiers $x_1 = \theta_1$ and $x_2 = \theta_2$ (Dmitruk and Kaganovich, 2008). Through its application, which is not fully detailed in this paper, we obtain that (OCP) admits an optimal control which can be defined as a very simple feedback.

**Theorem 1.** The optimal strategy $u(x)$ solution of (OCP) for B-admissible trajectories is the feedback control

$$ u(x) = \begin{cases} u_{min} & \text{if } x \in (S_{u_{max}})^-, \\ u_{max} & \text{if } x \in (S_{u_{max}})^+ \cup (S_{u_{max}}). \end{cases} $$

Note that $u(x)$ is not defined in $B_{11}$ due to the lack of control in the region. Additionally, for the case starting in $(S_{u_{max}})^-$, one has the following result.

**Proposition 2.** For initial conditions in $(S_{u_{max}})^-$,

* the optimal control consists of two bang arcs, as the suboptimal control (2);

* the optimal trajectories passes by the unstable Filippov equilibrium $(\theta_1, \theta_2)$.

The proof of such result involves showing there are no singular arcs in the optimal control, and thus $u(t)$ can only be a concatenation of bang arcs. Additionally, because of the two-dimensional affine structure in each regular domain, the sign of the switching function in the Hamiltonian can switch at most once throughout the whole interval $[0, t_f]$. Consequently, the optimal control consists of at most two bang arcs ($u_{min}$ or $u_{max}$), and the problem is reduced to finding the optimal switching time between the two arcs.
## 4.1 Numerical results

We illustrate our results with numerical simulations performed with Bocop (Team Commands, 2017), an open-source toolbox for solving OCPs<fnref n="3" />. The original problem is solved through a direct method, by approximating it by a finite dimensional optimization problem, using a Lobato time discretization method. The algorithm requires $s^-$ to be regularized to a smooth function. We then define, for $x \in \mathbb{R}$ and $k \in \mathbb{N}$, the Hill function

$$ \delta(x_i, \theta_i, k) = \frac{\theta_i^k}{x_i^k + \theta_i^k}, \quad (3) $$

which can approximate $s^-$ for large values of $k$ and, when $k \to \infty$,

$$ \lim_{k \to \infty} \delta(x_i, \theta_i, k) = \begin{cases} 1 & x_i < \theta_i, \\ 0 & x_i > \theta_i, \\ 1/2 & x_i = \theta_i. \end{cases} $$

Replacing $s^-$ by Hill functions (3) in system (S) yields the non-hybrid system

$$ \begin{cases} \dot{x}_1 = -\gamma_1 x_1 + u k_1 \delta(x_2, \theta_2, k), \\ \dot{x}_2 = -\gamma_2 x_2 + u k_2 \delta(x_1, \theta_1, k). \end{cases} $$

System parameters are fixed to $\gamma_1 = 1.2$, $\gamma_2 = 2$, $\theta_1 = 0.6$, $\theta_2 = 0.4$ and $k_1 = k_2 = 1$, which verify Assumption 1; and control bounds are set to $u_{min} = 0.5$ and $u_{max} = 1.5$ satisfying Assumption 2. The parameter $k$ of the Hill function is set to $k = 500$, which proved an acceptable approximation of the $s^-$ function. Figure 2 shows an optimal trajectory representing the transition (high, low) to (low, high). As predicted in Proposition 2, the optimal control is a bang-bang control: it consists of a first phase $[0, t_s]$ of low synthesis control $u_{min}$ until $x$ reaches the separatrix $(S_{u_{max}})$, followed by a phase $[t_2, t_f]$ of high synthesis control $u_{max}$ until $x_2$ reaches $x_2^f$. As it is customary when solving OCPs with direct methods, the algorithm does not count on any *a priori* information of the structure of the optimal control. Yet, the obtained trajectory is in agreement with Theorem 1, which confirms our theoretical results. Moreover, the solver is not restricted to consider only B-admissible trajectories, which suggests that the solution found in this work is optimal not only for Problem (OCP) along B-admissible trajectories but also for the general (OCP), without imposing the domain sequences.

Figure 3 shows different trajectories starting from $(S_{u_{max}})^+$ and $(S_{u_{max}})^-$. The streamplot represents the closed-loop dynamics for the optimal control defined in Theorem 1. All trajectories starting in $(S_{u_{max}})^-$ approach asymptotically the point $\Phi^*$ (denoted by a cross) until they reach the separatrix, point at which the state slides over it towards the Filippov equilibrium $(\theta_1, \theta_2)$. The optimal control for trajectories starting in $(S_{u_{max}})^+$ consists in $u \equiv u_{max}$ for the whole interval $[0, t_f]$, and do not pass by the Filippov equilibrium.

**Remark 1.** As already mentioned in the introduction, the dynamics is not uniquely defined at the undifferentiated point $(\theta_1, \theta_2)$, and the proposed solution is obtained by making a choice of dynamics at this point. Hence, con-

<figure><img src="image_2.png" /><figcaption>Fig. 2. Optimal trajectory with $x_1^0 = 0.8$, $x_2^0 = 0.3$ and $x_2^f = 0.7$. Times $t_1$ and $t_2$ are the transition times at which the state meets $x_1(t_1) = \theta_1$ and $x(t_2) = (\theta_1, \theta_2)$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 3. Optimal trajectories starting from different initial points, with $x_2^f = 0.7$. The streamplot represents the vector field resulting from applying the optimal bang-bang strategy from Theorem 1.</figcaption></figure>

cerning a biological implementation of our time-optimal strategy, it seems more reasonable to apply $u(t) \equiv u_{min}$ during a slightly longer time $\dot{t}_s = t_s + \epsilon$ with a small $\epsilon > 0$.

## 5. CONCLUSION

We showed novel results addressing the time-optimal control problem of a bistable gene-regulatory network. Our work is inspired by the genetic toggle switch, a milestone in Systems Biology, and of great interest for engineering genetic circuits. While in previous works (Chaves and Gouzé, 2011), the bang-bang nature of the control is imposed as a constraint, we showed that such characteristic is indeed a requirement in order to produce minimum-time transitions. Results also indicate that optimal trajectories should pass by Filippov equilibrium corresponding to the undifferentiated state $(\theta_1, \theta_2)$, which, as already signaled, represents a highly relevant state from the biological point of view. We expect that our result could be generalized to more complex genetic regulatory networks, where it often occurs that trajectories belonging to a given domain may bifurcate in different domains, similarly to what happens in the toggle switch case.

<footnote n="3">³ In order to guarantee the reproducibility of the numerical results, the computations can be re-runned through an online version of Bocop available at https://ct.gitlabpages.inria.fr/gallery/bistable/bistable.html, as part of the control toolbox project.</footnote>
ACKNOWLEDGEMENTS

We are very grateful to Jean-Baptiste Caillau (Université Côte d'Azur / Inria, McTao team), Jean-Luc Gouzé (Inria, Biocore team) and Madalena Chaves (Inria, Biocore team) for enlightening discussions on the topic.

REFERENCES

Balázsi, G., van Oudenaarden, A., and Collins, J.J. (2011). Cellular decision making and biological noise: from microbes to mammals. *Cell*, 144(6), 910–925.

Boscain, U., Grönberg, F., Long, R., and Rabitz, H. (2014). Minimal time trajectories for two-level quantum systems with two bounded controls. *Journal of Mathematical Physics*, 55(6), 062106.

Boscain, U. and Mason, P. (2006). Time minimal trajectories for a spin 1/2 particle in a magnetic field. *Journal of Mathematical Physics*, 47(6), 062101.

Chambon, L., Belgacem, I., and Gouzé, J.L. (2020). Qualitative control of undesired oscillations in a genetic negative feedback loop with uncertain measurements. *Automatica*, 112, 108642.

Chaves, M. and Gouzé, J.L. (2011). Exact control of genetic networks in a qualitative framework: the bistable switch example. *Automatica J. IFAC*, 47(6), 1105–1112.

Dmitruk, A.V. and Kaganovich, A.M. (2008). The hybrid maximum principle is a consequence of pontryagin maximum principle. *Systems & Control Letters*, 57(11), 964–970.

Farcot, E. and Gouzé, J.L. (2008). A mathematical framework for the control of piecewise-affine models of gene networks. *Automatica J. IFAC*, 44(9), 2326–2332.

Filippov, A.F. (1960). Differential equations with discontinuous right-hand side. *Matematicheskii sbornik*, 93(1), 99–128.

Gardner, T., Cantor, C., and Collins, J. (2000). Construction of a genetic toggle switch in escherichia coli. *Nature*, 403, 339–42. doi:10.1038/35002131.

Glass, L. and Pasternack, J.S. (1978). Prediction of limit cycles in mathematical models of biological oscillations. *Bull. Math. Biology*, 40(1), 27–44. doi: 10.1007/bf02463128. URL https://doi.org/10.1007/bf02463128. Papers presented at the Society for Mathematical Biology Meeting (Univ. Pennsylvania, Philadelphia, Pa., 1976).

Lugagne, J.B., Carrillo, S.S., Kirch, M., Köhler, A., Batt, G., and Hersen, P. (2017). Balancing a genetic toggle switch by real-time feedback control and periodic forcing. *Nature communications*, 8(1), 1–8.

Mallet-Paret, J. and Smith, H.L. (1990). The poincaré-bendixson theorem for monotone cyclic feedback systems. *Journal of Dynamics and Differential Equations*, 2(4), 367–421.

Team Commands, I.S. (2017). Bocop: an open source toolbox for optimal control. http://bocop.org.

Trélat, E. and Zuazua, E. (2015). The turnpike property in finite-dimensional nonlinear optimal control. *Journal of Differential Equations*, 258(1), 81–114.

Appendix A. HYBRID MAXIMUM PRINCIPLE FOR TIME OPTIMAL CONTROL

We adapt here the version of the Hybrid Maximum Principle given in Dmitruk and Kaganovich (2008) to the time optimal setting.

Let $t_0 < t_1 < \dots < t_\nu$ be real numbers. Denote by $\Delta_k$ the time interval $[t_{k-1}, t_k]$. For continuous functions $x^k: [t_0, t_\nu] \to \mathbb{R}^n$, $k \in \{1, \dots, \nu\}$, define the vector
$$p = (t_0, (t_1, x^1(t_0), x^1(t_1)), \dots, (t_\nu, x^\nu(t_{\nu-1}), x^\nu(t_\nu))) \in \mathbb{R}^d,$$
where $d = 1 + (2n + 1)\nu$. Let $(f_k)_{k \in \{1, \dots, \nu\}}$ be smooth vector fields on $\mathbb{R}^n$, and $(\phi_i)_{i \in \{1, \dots, m\}}$, $(\eta_j)_{j \in \{1, \dots, q\}}$ be two families of smooth functions defined on $\mathbb{R}^{(\nu+1)(n+1)}$. For $t \in [t_0, t_\nu]$ and a collection $(U_k)_{k \in \{1, \dots, \nu\}}$ of subsets of $\mathbb{R}^q$, $q \ge 1$, consider the autonomous hybrid optimal control problem

$$
\left\{
\begin{array}{l}
\displaystyle \text{minimize } t_{\nu} - t_0 \\
\dot{x}^k(t) = f_k(x^k(t), u^k(t)), \\
u^k(t) \in U_k, \\
t \in \Delta_k, k \in \{1, \dots, \nu\} \\
\eta_j(p) = 0, j = 1, \dots, q \\
\phi_i(p) \le 0, i = 1, \dots, m.
\end{array}
\right.
\qquad \text{(HYBRID)}
$$

Definition 2. For a tuple $w = (t_0; t_k, x^k(t), u^k(t), k = 1, \dots, \nu)$ which is extremal for Problem (HYBRID), define:

• the trajectory $(x(t))_{t \in [t_0, t_\nu]}$ which is equal to $x^k(t)$ for every $t \in \Delta_k \setminus \{t_k\}$ and $k \in \{1, \dots, \nu\}$;

• the adjoint trajectory $(\lambda(t))_{t \in [t_0, t_\nu]}$ which is equal to $\lambda^k(t)$ for every $t \in \Delta_k \setminus \{t_k\}$ and $k \in \{1, \dots, \nu\}$;

• the control $(u(t))_{t \in [t_0, t_\nu]}$ which is equal to $u^k(t)$ for every $t \in \Delta_k \setminus \{t_{k-1}\}$ and $k \in \{1, \dots, \nu\}$.

Define, for every $k \in \{1, \dots, \nu\}$, $t \in \Delta_k$,

$$
H^k(x^k, \lambda^k, \lambda_0, u^k) = \langle \lambda^k, f_k(x^k, u^k) \rangle - \lambda_0.
$$

Theorem 3. Assume that $(\tilde{x}(\cdot), \tilde{u}(\cdot), \tilde{p})$ is an optimal solution of Problem (HYBRID). Then there exists

$$
(\alpha, \beta, \lambda(\cdot), \lambda_0)
$$

, where $\alpha = (\alpha_1, \dots, \alpha_m) \in \mathbb{R}^m$, $\beta = (\beta_1, \dots, \beta_q) \in \mathbb{R}^q$, $\lambda = (\lambda^1, \dots, \lambda^\nu)$, all $\lambda^k : \Delta_k \to \mathbb{R}^n$ for $k \in \{1, \dots, \nu\}$ being Lipschitz functions, and a constant $\lambda_0 \ge 0$ such that:

• $(\lambda_0, \alpha, \beta) \neq 0$;

• For every $i \in \{1, \dots, m\}$, $\alpha_i \ge 0$;

• For every $i \in \{1, \dots, m\}$, $\alpha_i \phi_i(\tilde{p}) = 0$;

• For almost every $t \in \Delta_k$,

$$
\begin{align*}
\dot{x}^k &= \frac{\partial H^k}{\partial \lambda}(x^k, \lambda^k, \lambda_0, \tilde{u}) \\
\dot{\lambda}^k &= -\frac{\partial H^k}{\partial x^k}(x^k, \lambda^k, \lambda_0, \tilde{u}) \tag{E} \\
H^k(x^k, \lambda^k, \lambda_0, \tilde{u}) &= \max_{u \in \Omega} H^k(x^k, \lambda^k, \lambda_0, u) = 0.
\end{align*}
$$

Moreover, if we define $L(p) = \lambda_0(t_\nu - t_0) + \sum_{i=1}^{m} \alpha_i \phi_i(p) + \sum_{j=1}^{q} \beta_j \eta_j(p)$, then we have the following transversality and discontinuity conditions at times $t = t_0, \dots, t_\nu$:

• At the initial and final times $t_0$ and $t_\nu$, we have
$$ \lambda^1(t_0) = \frac{\partial L}{\partial x^1(t_0)}(\tilde{p}) $$

$$ \lambda^{\nu}(t_{\nu}) = \frac{\partial L}{\partial x^{\nu}(t_{\nu})}(\tilde{p}) $$

* At the crossing times $(t_k)_{k \in \{1, \dots, \nu-1\}}$, we have, for every $k \in \{1, \dots, \nu-1\}$,

$$ \lambda^k(t_{k-1}) = \frac{\partial L}{\partial x^k(t_{k-1})}(\tilde{p}) $$

$$ \lambda^k(t_k) = -\frac{\partial L}{\partial x^k(t_k)}(\tilde{p}). $$