# On $L^1$ and time-optimal state transitions in piecewise linear models of gene-regulatory networks

Agustín G. Yabo, Nicolas Augier

<figure><img src="image_1.png" /></figure>

► To cite this version:

Agustín G. Yabo, Nicolas Augier. On $L^1$ and time-optimal state transitions in piecewise linear models of gene-regulatory networks. Nonlinear Analysis: Hybrid Systems, 2026, 59, pp.101650. 10.1016/j.nahs.2025.101650. hal-04820387v3

HAL Id: hal-04820387

https://hal.science/hal-04820387v3

Submitted on 14 Nov 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

ELSEVIER

Contents lists available at ScienceDirect

Nonlinear Analysis: Hybrid Systems

journal homepage: www.elsevier.com/locate/nahs

<figure><img src="image_4.png" /></figure>

# On $L^1$ and time-optimal state transitions in piecewise linear models of gene-regulatory networks☆

<figure><img src="image_5.png" /></figure>

Check for
updates

Agustín G. Yabo a,*, Nicolas Augier b,1

## ARTICLE INFO

### Keywords:

- Hybrid optimal control

- Mathematical biology

- Gene regulatory networks

- Piecewise linear systems

- Hybrid systems

## ABSTRACT

In this paper, we investigate optimal state transfers for a generic class of piecewise-linear models widely used to qualitatively describe gene-regulatory networks. Motivated by the main practical drawbacks of artificially regulating gene expression through chemical inducers, the optimality of the transitions is defined as the convex combination of the total time and the $L^1$ cost of the control. Solutions are studied through a Hybrid Pontryagin's Maximum Principle approach, which allows to characterize the optimal trajectories and control for the general formulation of the problem. Then, we focus on two practical examples of two-dimensional regulatory networks: the bistable switch, for which the objective is to induce optimal transitions between its two stable steady states, and the damped genetic oscillator, where the goal is to induce sustained oscillatory behaviors. The resulting optimal control strategies can be expressed in state feedback form, involving both bang arcs and inactive control periods, and are shown to slide over certain separatrices of the uncontrolled system that characterize the boundaries of the admissibility set.

## 1. Introduction

Most of the underlying biological processes in living beings are governed by the numerous biochemical interactions between genes occurring at the intracellular level. Even for the simplest natural phenomena, these interactions can form complex regulatory networks exhibiting highly nonlinear behaviors. For decades, the study of these mechanisms, often referred to in the literature as Gene Regulatory Networks (GRNs), has been crucial in understanding elemental cellular functions such as signaling, differentiation, metabolism and growth. In the context of systems biology, a wide spectrum of quantitative and qualitative modeling paradigms have been used to describe and study GRNs, ranging from logic models to stochastic differential equations [1]. These approaches are often used in a complementary manner, as they can respond to different research questions, and thus provide insight into different biological mechanisms. In systems theory, it is customary to describe the expression levels of genes and proteins as solutions of deterministic differential equations, and to model their interactions as nonlinear mathematical expressions of chemical kinetics. A particular feature of GRNs is that they operate in conditions in which synthesis rates are almost always saturated [2]. Moreover, and in spite of the simplicity of such models, the nonlinearities can often obscure the study of a given phenomenon, even if they do not play a decisive role in the behavior of the system. This motivated the development of the piecewise linear (PWL) framework [3,4], a particular hybrid/nonlinear modeling paradigm in which the nonlinear gene expression rates are approximated through piecewise constant functions.

<footnote>ª MISTEA, Université Montpellier, INRAE, Institut Agro, Montpellier, France</footnote>

<footnote>ᵇ Mac Team, LAAS-CNRS, Toulouse, France</footnote>

<footnote>☆ This article is part of a Special issue entitled: 'Nonsmooth Dynamical Systems' published in Nonlinear Analysis: Hybrid Systems.</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: agustin.yabo@inrae.fr (A.G. Yabo).</footnote>

<footnote n="1">¹ These authors contributed equally to this work. The order of author names was randomly determined.</footnote>
The PWL modeling framework has been instrumental in characterizing the qualitative behavior of numerous naturally-evolved regulatory processes. For instance, it has been used to study the transitions between two genetic states (normal vs. starvation) in the bacterium *Bacillus subtilis* [5]. Furthermore, this perspective has raised interesting theoretical questions for the control systems community, in particular in how to modify these natural behaviors through an external action to yield specific dynamical patterns [6,7]. Similarly, more recent works have also studied how to synchronize networks of coupled GRNs [8,9]. The aforementioned approaches focus mostly on proposing feasible control laws, and so they do not consider the optimality of these control strategies, which is an aspect that has received significantly less attention from the community. For instance, the time needed to induce transfers between states can be a major constraint in certain biosynthetic devices that involve multiple asynchronous mechanisms functioning at different timescales [10]. In this sense, time-efficiency of state transfers in GRNs becomes essential for biotechnological applications, thus motivating an optimal control perspective [11,12].

Optimality of control strategies can be defined in regard to numerous criteria, depending on the chosen control scheme and objectives. In an experimental setting, gene networks can be externally regulated by activating or inactivating gene expression through chemical inducibles. The most classical example is the use of the diffusible molecule IPTG (isopropyl β-d-thiogalactoside) to induce protein expression in E. coli, a method widely used, among others, to externally control the state of biosynthetic devices [13,14]. IPTG acts on the *lac* operon—which is a group of genes used by bacteria to metabolize lactose—by inhibiting the gene *lacI* that blocks gene transcription in the absence of lactose (which consequently triggers the transcription of the *lac* operon). The main drawback of modulating gene expression with this method is that such chemical inducers are highly expensive [15] and produce toxic elements known to induce stress to cells [16]. In this context, it is compulsory to limit the control action when designing efficient feedback strategies in GRNs, which can be done either by imposing upper bounds onto the total control usage, or by considering it as an additional operation cost.

This is the question addressed in this work, where we study optimal state transfers for a generic class of PWL mathematical models of GRNs, in terms of the cost associated to the use of control and the duration of the transition. For that, we consider a convex combination of the $L^1$ norm of the control function and the total transfer time. The mathematical originality of the approach lies in the integration of a hybrid framework on the state with a non-smooth mixed $L^1$ and time-OCP (Optimal Control Problem). On one hand, the exploration of $L^1$ cost functions has emerged as a significant research area [17–19] characterized by the presence of inactivated arcs (i.e. arcs along which the cost is zero) in optimal synthesis, as a generic feature of optimal control solutions [20]. Substantial efforts have also been directed towards understanding the optimality of singular arcs in this particular setting [17]. On the other hand, the optimal control of hybrid systems has increasingly received attention in recent years, with comprehensive works addressing optimality conditions for this particular framework [21,22]. However, very few studies have applied these results to practical problems from applied sciences, let alone mathematical and synthetic biology. Recent contributions [23,24] aim to clarify the application of the Hybrid Maximum Principle (HMP) across various frameworks, striving to make it a readily usable tool for non-smooth control systems defined regionally. A difficulty is that complex behaviors can occur at the interfaces where dynamic switches happen between different regions. In this work, we mitigate the latter behaviors by employing a sequential-in-time HMP as introduced in [25], effectively disregarding certain trajectory crossings at highly singular interfaces defined by the thresholds of the gene network under consideration.

In addition to the results obtained for the general formulation, we investigate the problem of inducing efficient state transitions in two different examples of GRNs, both widely studied in the field of systems biology. The latter systems represent the smallest GRNs exhibiting bistability and oscillations, and can be modeled through two-dimensional positive and negative feedback loops, respectively. The first case is a bistable system inspired on the genetic toggle switch, a biosynthetic flip-flop device consisting of two genes mutually repressing each other, for which the control objective is to induce a switch between states. The genetic toggle switch was first implemented in the bacterium E. coli by exploiting the mutual repression between the genes *lacI* and *tetR*, whose concentrations account for the state of the switch that can be controlled by externally modulating their synthesis rates through the inducers IPTG and aTc [26]. This device represents a milestone in biocomputing, for being the first implementation of a biosynthetic memory unit capable of storing 1 bit of information in a living cell. The second example is a genetic damped oscillator, for which the control objective is to efficiently transform the damped oscillations into sustained ones. Genetic oscillators have been extensively used in the literature to qualitatively describe the dynamical behavior of more complex intracellular phenomena [27], such as the mammalian circadian clock, or the carbon starvation response in E. coli [28]. For each case, the biological objective can be represented through a state transition control problem, for which we investigate the structure of the optimal solutions by applying the developed theoretical framework and, in parallel, by exploiting the inherent dynamical features of each individual system, and the reachability properties in terms of the control and the chosen terminal state. A thorough investigation of the problems shows that optimal state transfers can be expressed as state-feedback laws that, when applied to the systems, yield piecewise-constant control functions involving both bang (i.e. saturated) and inactive control periods. Optimal trajectories are shown to slide over certain surfaces of the uncontrolled system (that we call separatrices), which are specific to the dynamics of each system, and characterize the boundaries of the admissibility set. Moreover, we show how the interplay between time and $L^1$-optimality reveals the existence of bifurcations among the different possible control strategies, depending on the weight assigned to each objective.

The paper is organized as follows: In Section 2, we define the class of systems of interest and the optimal control problem, followed by the general results concerning the structure of the optimal synthesis, through the use of an adaptation of the Hybrid Maximum Principle stated in Theorem 1. Then, in Section 3, we focus on the two example applications, where the general theoretical results are refined to each particular case. For each problem, an optimal feedback control law is derived, and the results are verified through numerical simulations. Note that the theoretical results of Section 2 are independent of the interaction graph of each system, and could be used for the understanding of more complex motifs.
<figure><img src="image_6.png" /><figcaption>Fig. 1. Streamplot of a bidimensional non-controlled (i.e. $F_2(x) = 0$) sample system with multiple thresholds $\theta$, illustrating the PWL dynamics that can be obtained in such systems, and the resulting multistability patterns.</figcaption></figure>

## 2. General framework

In this section, we define a general mathematical model of an n-dimensional GRN. For $i \in \{1, 2\}$ and $j \in \{1, ..., n\}$, let $f_j^i$ be given functions defined on $\{0, 1\}^n$ representing the (catalytic or inhibitory) influence of other genes onto gene j, and define the sequences $(\theta_{q,j}^i)_{q,j \in \{1,...,n\}}$ of non-negative numbers such that for given $q, j \in \{1, ..., n\}$, the value $\theta_{q,j}^i$ corresponds to the transcriptional threshold for the influence of gene j on gene q. Define for every $x \in \mathbb{R}^n$, $\kappa_q^i(x) = (s^\pm(x_1, \theta_{q,1}^i), ..., s^\pm(x_n, \theta_{q,n}^i))$, and consider the PWL dynamical model denoted by

$$ \dot{x} = F_1(x) + u(t)F_2(x) - \Gamma x \quad (1) $$

where $F_1$ and $F_2$ are the vectors of (non-controlled and controlled, respectively) synthesis rates, defined as

$$ F_1(x) = \begin{bmatrix} f_1^i(\kappa_1^i(x)) \\ \vdots \\ f_n^i(\kappa_n^i(x)) \end{bmatrix}, $$

and $u(t) \in [u_{\min}, u_{\max}]$ is the external control function, with bounds satisfying $0 < u_{\min} < 1 < u_{\max}$. The state $(x(t))_{t \ge 0}$ represents the concentrations of each protein, while the matrix $\Gamma$ is positive diagonal with pairwise distinct coefficients $(\gamma_j)_{j \in \{1,...,n\}}$, and for $j \in \{1, ..., n\}$, $y \in \mathbb{R}$, and for $\theta \in \mathbb{R}$, the functions $s^+(\cdot, \theta) : \mathbb{R} \to \mathbb{R}$ and $s^-(\cdot, \theta) : \mathbb{R} \to \mathbb{R}$ are such that

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

For simplicity of notations, we denote by $s^\pm$ a function which is either $s^+$ or $s^-$. It is assumed that $s^\pm(y, \theta) \in [0, 1]$ for $y = \theta$. For a fixed $j \in \{1, ..., n\}$, some different thresholds may be involved in the dynamics for the same variable $x_j$. Next, we state some dynamical properties of the general system; we refer to Fig. 1 for an example of such dynamical behavior. Since $(x(t))_{t \ge 0}$ represents protein concentrations, we formulate the following assumption:

**Assumption 1.** For every constant input $u(t) \equiv u \in [u_{\min}, u_{\max}]$, the quadrant $(\mathbb{R}_+)^n$ is forward invariant by the dynamics of Eq. (1).

**Definition 1.** Define the regular domains as non-empty finite intersections over $j, q \in \{1, ..., n\}$ of domains of the type

$$ B_{j,q}^{i,+} = \{x \in \mathbb{R}^n \mid x_j > \theta_{j,q}^i\}, $$

or

$$ B_{j,q}^{i,-} = \{x \in \mathbb{R}^n \mid x_j < \theta_{j,q}^i\}, $$

for $i \in \{1, 2\}$, in which the vector fields $F_1$ and $F_2$ are constant. Note that an open subset $B$ of $\mathbb{R}^n$ is a regular domain if and only if it can be written as a non-empty set under the form

$$ B = \bigcap_{i \in \{1, 2\}} \bigcap_{j,q \in \{1, ..., n\}} B_{j,q}^{i, \epsilon_{i,j,q}}, $$

where $\epsilon_{i,j,q} \in \{+, -\}$ for every $i \in \{1, 2\}$ and $j, q \in \{1, ..., n\}$.

**Definition 2.** Eq. (1) restricted to a regular domain $B$ is a linear dynamical system on $\mathbb{R}^n$ having an asymptotically stable equilibrium $\phi$, called focal point. If $\phi \in B$, then $\phi$ is a steady state for Eq. (1).
The control $u$ which is scalar and real valued, is assumed to belong to the functional space $L^\infty(\mathbb{R}, [u_{\min}, u_{\max}])$, so that we can define the solutions of Eq. (1) in Filippov sense. In a given admissible domain $B$, the vector fields $F_1(x)$ and $F_2(x)$ are constant vectors of $\mathbb{R}^n$.

## 2.1. Optimal control problem

We consider the problem of inducing a state transfer for system (1) under the general initial and terminal conditions:

$$x(0) = x_0 > 0, \quad x(t_f) \in \mathcal{K} \subset B_f,$$

where $x_0$ belongs to a regular domain, and $\mathcal{K}$ is the restriction of an affine semialgebraic set to the final regular domain $B_f$. For a fixed $\lambda \in (0, 1]$, the cost function to minimize is defined as

$$J_{\lambda}(u) = \lambda \int_{0}^{t_f} |u(t) - 1| dt + (1 - \lambda)t_f$$

among admissible trajectories of Eq. (1), and $t_f \ge 0$ which is a free non-negative real number. In order to solve this optimal control problem, one has to face two difficulties:

* The hybrid nature of the dynamics;

* The non-smoothness of the $L^1$ Lagrangian cost.

The first point can be tackled by optimizing over $B$-admissible trajectories, as defined in the next definition, while the second point can be solved by applying the Hybrid Pontryagin's Principle from [25] to an extended system, as described in Section 2.2.

**Definition 3.** Let $\nu$ be a non-negative integer, and a sequence $B = (B_k)_{k \in \{1, \dots, \nu\}}$ of regular domains as defined in **Definition 1**.

* We say that a solution $x(t)$ of Eq. (1) is $B$-admissible if it is absolutely continuous and there exists a time $t_f > 0$, a control $u(\cdot) \in L^\infty([0, t_f], [u_{\min}, u_{\max}])$, and times $t_0 = 0 < t_1 < \dots < t_\nu = t_f$ such that $x(t) \in B_k$ for every $t \in \Delta_k$, where $\Delta_k = (t_k, t_{k+1})$, and define $x^k : \Delta_k \to \mathbb{R}^n$ as the continuous function which is the restriction of $x(\cdot)$ to the time interval $\Delta_k$, and $u^k : \Delta_k \to \mathbb{R}$ as the $L^\infty$ function which is the restriction of $u(\cdot)$ to $\Delta_k$.

* Define the vector $z \in \mathbb{R}^d$ as

$$z = (t_0, (t_1, x^1(t_0), x^1(t_1)), \dots, (t_\nu, x^\nu(t_{\nu-1}), x^\nu(t_\nu))),$$

where $d = 1 + (2n + 1)\nu$, and parametrize the intersection wall between the closures $\bar{B}_j$ and $\bar{B}_{j+1}$ as a finite intersection of sets under the form $\{\phi_j^k(z) = 0\} \cap \{\eta_j^k(z) \ge 0\}$, for affine functions $(\phi_j^k, \eta_j^k)_{j \in \{1, \dots, q_k-1\}}$ depending on $z$, for every $k \in \{1, \dots, \nu-1\}$, where $q_k \ge 1$, and $t_0 = 0$ and $t_\nu = t_f$.

* The set $\mathcal{K}$ is parametrized as the intersection of an affine semialgebraic set with $B_f$, hence it can be written

$$\mathcal{K} = \bigcap_{j=1}^{q_\nu} \{\phi_j^\nu(z) = 0\} \cap \{\eta_j^\nu(z) \ge 0\},$$

for given affine functions $(\phi_j^\nu, \eta_j^\nu)_{j \in \{1, \dots, q_\nu\}}$, where $q_\nu \ge 1$.

Under this definition, the optimization is performed over the set of trajectories following a predetermined sequence of regular domains $B$, and the continuity of trajectories is in practice imposed by adding some supplementary constraints, as it will be illustrated in **Example 1**. Of course, this implies that the approach is restricted to biological systems—and any other system that can be represented through PWL models—whose dynamical behavior is known beforehand. This is normally the case when the objectives correspond to biologically-inspired problems, such as toggling a genetic switch or sustaining oscillations, which are the examples studied in Section 3. However, grasping the dynamical specificities of higher-dimensional problems can become a more difficult task, for which other tools can be useful (see the beginning of Section 3 for a discussion on dedicated numerical methods).

**Example 1.** In the simple two-dimensional case where all the synthesis rates are controlled, we have $n = 2$, $\nu = 2$, $F_1 \equiv 0$, and

$$F_2(x) = \begin{pmatrix} f_1(s^{\pm}(x_1, \theta_1), s^{\pm}(x_2, \theta_2)) \\ f_2(s^{\pm}(x_1, \theta_1), s^{\pm}(x_2, \theta_2)) \end{pmatrix},$$

with thresholds $\theta_1, \theta_2 > 0$. This defines the regular domains

$$B_1 = \{x \in \mathbb{R}^2 \mid x_j < \theta_j, j \in \{1, 2\}\},$$

$$B_2 = \{x \in \mathbb{R}^2 \mid x_1 < \theta_1, x_2 > \theta_2\},$$

and so we can define $B = (B_1, B_2)$-admissible trajectories starting in $x(t_0) = (x_1^0, x_2^0) \in B_1$ and ending in $x(t_f) = (x_1^f, x_2^f) \in B_2$ with the constraint functions defined, for $z = (t_0, (t_1, x^1(t_0), x^1(t_1)), (t_2, x^2(t_1), x^2(t_2))) \in \mathbb{R}^{11}$, as $t_0 = 0$, $x_1^1(t_0) - x_1^0 = 0$, $x_2^1(t_0) - x_2^0 = 0$, $t_2 - t_f = 0$, and

$$\eta_1^1(z) = \theta_1 - x_1^1(t_1) \ge 0,$$
$$\phi_1^1(z) = x_2^1(t_1) - \theta_2 = 0,$$

$$\phi_1^2(z) = x_2^2(t_1) - \theta_2 = 0,$$

$$\phi_1^3(z) = x_1^2(t_2) - x_1^f = 0,$$

$$\phi_2^3(z) = x_2^2(t_2) - x_2^f = 0,$$

$$c_1^2(z) = x_1^2(t_1) - x_1^1(t_1) = 0,$$

where the constraint $c_1^2(z) = 0$ is added in order to ensure the continuity of the trajectories, as required by Definition 3.

## 2.2. Adapted hybrid Pontryagin's Principle

Let $v \ge 1$, and a sequence $B = (B_j)_{j \in \{1, \dots, v\}}$ of regular domains as defined in Definition 1, and associated affine constraint functions $(\phi_j, \eta_j)_j$. In the notations of Definition 3, define, for every $k \in \{1, \dots, v\}$, $x^k \in B_k$, $u^k \in [u_{\min}, u_{\max}]$ and $p^k \in \mathbb{R}^n$, the pseudo-Hamiltonian

$$\mathcal{H}^k(x^k, p^k, p^0, u^k) = \langle p^k, F_k(x^k, u^k) \rangle + p^0(\lambda|u^k - 1| + 1 - \lambda),$$

where $F_k$ is the restriction of the vector field $F_1(x) + uF_2(x) - \Gamma x$ for $x \in B_k$. For a given B-admissible trajectory $(x(t), u(t))$, define, for every $t \ge 0$, $y(t) = \lambda \int_0^t |u(s) - 1| ds + (1 - \lambda)t$, and define $Z^k(t) = (x^k(t), y^k(t))$, so that the extended system in the variable $Z^k(t)$ satisfies

$$\dot{Z}^k = \begin{pmatrix} F_k(x^k, u^k) \\ \lambda|u^k - 1| + (1 - \lambda) \end{pmatrix},$$

and minimizing $J_\lambda$ is then equivalent to minimizing $y(t_f)$. It is then possible to apply the temporally Hybrid Pontryagin's Maximum Principle (HMP) from [25, Theorem 4] (we can also refer to [29] for a comparison of the different settings for Hybrid OCPs), for which the minimization is made among trajectories following a fixed sequence of controlled vector fields, then restricting the minimization to B-admissible trajectories. We obtain the following necessary optimality conditions for the System (1). In order to improve the readability of the paper, we will replace the constraint functions $(\phi_j^k)_{j,k}$ and $(\eta_j^k)_{j,k}$ indexed by pairs of integers $(j, k)$ as in Definition 3 by two sequences $(\phi_j)_{j \in \{1, \dots, m\}}$ and $(\eta_j)_{j \in \{1, \dots, q\}}$ where $m, q \ge 1$. These functions incorporate the intermediate constraints at the walls between the regular domains and the terminal constraint $x(t_f) \in \mathcal{K}$; as well as the constraints imposing the continuity of trajectories at the intersection walls between regular domains. Note that our framework excludes sliding motions along the boundaries of regular domains.

**Theorem 1.** Let $(x(t), u(t))$ be a B-admissible trajectory associated with a time sequence $\{t_0 = 0, \dots, t_v = t_f\}$ which is optimal for the minimization of $J_\lambda(u)$. Then there exists

$$(\alpha, \beta, p(\cdot), p^0),$$

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

* $(p^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(z) = 0$;

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

$$\dot{x}^k = \frac{\partial \mathcal{H}^k}{\partial p}(x^k, p^k, p^0, u),$$

$$\dot{p}^k = - \frac{\partial \mathcal{H}^k}{\partial x^k}(x^k, p^k, p^0, u), \qquad (E)$$

$$\mathcal{H}^k(x^k, p^k, p^0, u) = \max_{v \in [u_{\min}, u_{\max}]} \mathcal{H}^k(x^k, p^k, p^0, v) = 0. \qquad (V)$$

Moreover, if we define the Lagrangian

$$L(z) = p^0 \left( \lambda \int_{t_0}^{t_v} |u-1|dt + (1-\lambda)(t_v-t_0) \right) + \sum_{j=1}^{m} \alpha_j \phi_j(z) + \sum_{j=1}^{q} \beta_j \eta_j(z),$$

then we have the following transversality and discontinuity conditions at times $t \in \{t_0 = 0, \dots, t_v = t_f\}$:

* At the initial and final times $t_0 = 0$ and $t_v = t_f$, we have

$$p^1(0) = -\frac{\partial L}{\partial x^1(0)}(z),$$

$$p^v(t_f) = -\frac{\partial L}{\partial x^v(t_f)}(z).$$
• At the crossing times $(t_k)_{k \in \{1, \dots, v-1\}}$, we have, for every $k \in \{1, \dots, v-1\}$,

$$p^k(t_{k-1}) = \frac{\partial L}{\partial x^k(t_{k-1})}(z), \qquad (J)$$

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

Analogously, we say that an extremal $(x(t), p(t), p^0, u(t))$ solution of HMP is B-admissible if its associated solution $(x(t), u(t))$ is B-admissible.

**Remark 1.** In the continuous framework, a process $(x(t), p(t), p^0, u(t))$ is said to be *admissible* if the trajectory is bounded and satisfies the terminal constraints. In the framework here described, an additional requirement is to comply with the B-admissibility condition, which allows the study of a space-stratified OCP through its associated time-partitioned problem as done in this section (see e.g. [29] for more details).

### 2.2.1. Adjoint equation and some useful quantities

For $x$ belonging to a given regular domain, $u \in [u_{\min}, u_{\max}]$, $p \in \mathbb{R}^n$ and $p^0 \le 0$, define the functions

$$\phi^{-}(x, p, p^{0}) = \langle F_{2}(x), p \rangle + \lambda p^{0},$$

$$\phi^{+}(x, p, p^{0}) = \langle F_{2}(x), p \rangle - \lambda p^{0},$$

that satisfy $\phi^+ \ge \phi^-$. With this notation, the Hamiltonian $\mathcal{H}^k$ in each regular domain becomes

$$\mathcal{H}^k(x, p, p^0, u) = \begin{cases} \langle F_1(x) - \Gamma x, p \rangle + u \phi^+(x, p, p^0) + p^0(1 - 2\lambda) & \text{if } u \ge 1, \\ \langle F_1(x) - \Gamma x, p \rangle + u \phi^-(x, p, p^0) + p^0 & \text{if } u < 1. \end{cases}$$

Applying Theorem 1, we obtain that along an extremal $(x(t), p(t), p^0, u(t))_{t \in [0, t_f]}$, the dynamics of the adjoint state (E) is

$$\dot{p}(t) = \Gamma p(t), \qquad (2)$$

for every $t \in [0, t_f]$, and it follows that

$$(\phi^{-})^{(m)}(t) = (\phi^{+})^{(m)}(t) = \langle F_{2}(x(t)), \Gamma^{m} p(t) \rangle,$$

for a.e. $t \in [0, t_f]$, and $m \ge 1$.

**Definition 4.** We say that a B-admissible extremal $(x(t), p(t), p^0, u(t))$ is normal (respectively, abnormal) if $p^0 \ne 0$ (respectively, $p^0 = 0$).

In accordance with the definition given in [17], we give the following definition of singular arcs.

**Definition 5.** We say that an arc $(x(t), p(t), p^0, u(t))_{t \in [t_1, t_2]}$ of a B-admissible extremal $(x(t), p(t), p^0, u(t))_{t \in [0, t_f]}$ is singular if it satisfies

$$\phi^{-}(x(t), u(t), p(t), p^{0})\phi^{+}(x(t), u(t), p(t), p^{0}) \equiv 0,$$

for a.e. $t \in [t_1, t_2]$.

A particular property of the state transfer problem for systems under the form (1) will be shown in Theorem 2, stating in particular that singular arcs cannot appear along normal trajectories.

### 2.3. General qualitative structure of optimal trajectories

Let $(B_k)_{k \in \{1, \dots, v\}}$ be a sequence of regular domains, and for $x \in B_k$, let $q = \#\Xi_k$, where $\Xi_k = \{j \in \{1, \dots, n\} \mid \langle F_2(x), e_j \rangle \ne 0\}$, and $\#\Xi_k$ denotes the cardinality of $\Xi_k$, and let $K_2^k \in \mathbb{R}^q$ be the vector of non-zero components of $F_2(x)$. The following result states that every normal optimal trajectory is a concatenation of bang arcs (i.e. $u(t) \in \{u_{\min}, u_{\max}\}$) and inactivated arcs (i.e. $u(t) \equiv 1$). Its proof can be found in Appendix B.

**Theorem 2.** Along a normal extremal, the optimal control $u(t)$ is non-singular, and it takes values in the set $\{u_{\min}, 1, u_{\max}\}$ for a.e. $t \in [0, t_f]$.

**Remark 2.** For the minimal-time problem (i.e., $\lambda = 0$), there are in general no inactivated arcs along extremal trajectories excepted when they admit singular arcs which might coincide with inactivated arcs.
<figure><img src="image_7.png" /><figcaption>Fig. 2. Examples of functions $\phi^{\pm}(t)$ along normal (i.e. $p^0 \neq 0$) and abnormal (i.e. $p^0 = 0$) trajectories with no singular arcs.</figcaption></figure>

## 2.4. Bounds on the number of switches for normal trajectories

Assume in this section that $\lambda \in (0, 1]$ and $p^0 \neq 0$, so that the conclusion of Theorem 2 holds.

**Proposition 3.** Let $(x(t), p(t), p^0, u(t))$ be a normal extremal. Then, for given $t_2 > t_1 \ge 0$, along an arc in the subinterval $[t_1, t_2]$, the optimal control is

$$u(t) = \begin{cases} u_{\max} & \text{if } \phi^{-}(x(t), p(t), p^{0}) > 0, \\ u_{\min} & \text{if } \phi^{+}(x(t), p(t), p^{0}) < 0, \\ 1 & \text{if } \phi^{+}(x(t), p(t), p^{0}) \ge 0 \text{ and } \phi^{-}(x(t), p(t), p^{0}) \le 0. \end{cases}$$

**Proof.** According to the maximization condition of Theorem 1, the control $u$ should maximize the Hamiltonian. If $\phi^{-}(x(t), p(t), p^0) > 0$ for $t \in [t_1, t_2]$, then the Hamiltonian $H^k(x(t), p(t), p^0, u)$ is maximized with $u(t) \equiv u_{\max}$ for a.e. $t \in [t_1, t_2]$. Analogously, if $\phi^{+}(x(t), u(t), p(t), p^0) < 0$ for $t \in [t_1, t_2]$, then the Hamiltonian $H^k(x(t), p(t), p^0, u_{\min})$ is maximized with $u(t) \equiv u_{\min}$ for a.e. $t \in [t_1, t_2]$. Finally, if $u(t) \equiv u_{\max}$ for $t \in [t_1, t_2]$ (respectively, $u(t) \equiv u_{\min}$), then $\phi^{-}(x(t), u(t), p(t), p^0) \ge 0$ (respectively, $\phi^{+}(x(t), u(t), p(t), p^0) \le 0$). Hence, arcs $u \equiv 1$ occur if and only if $\phi^{+}(x(t), u(t), p(t), p^0) \ge 0$ and $\phi^{-}(x(t), u(t), p(t), p^0) \le 0$. $\square$

In particular, in normal extremals there cannot be a direct switch between $u_{\min}$ and $u_{\max}$. We state the following classical property, which can be deduced by induction on $n \ge 1$ via a simple application of Rolle's Theorem.

**Lemma 1.** Let $(a_j)_{j \in \{0, \dots, n\}}$ be real numbers and $(\lambda_j)_{j \in \{1, \dots, n\}}$ be pairwise distinct real numbers. Then the function $\phi_n$ defined by $\phi_n(t) = a_0 + \sum_{j=1}^n a_j e^{\lambda_j t}$ for $t \in \mathbb{R}$ has at most $n$ zeros.

As a byproduct of Proposition 3 together with Lemma 1, we obtain the following result.

**Proposition 4.** Let B be a regular domain in the sense of Definition 1, and let $(x(t), u(t), p(t), p^0)$ be an extremal with associated trajectory satisfying $x(t) \in B$ for a.e. $t \in [t_1, t_2]$. Then the functions $t \mapsto \phi^{-}(x(t), u(t), p(t), p^0)$ and $t \mapsto \phi^{+}(x(t), u(t), p(t), p^0)$ vanish at most $n$ times in the interval $t \in [t_1, t_2]$.

**Corollary 1.** A normal optimal trajectory is made in each regular domains of:

* at most $n$ switches between $u(t) \equiv 1$ and $u(t) \equiv u_{\min}$ (independently of the order).

* at most $n$ switches between $u(t) \equiv 1$ and $u(t) \equiv u_{\max}$ (independently of the order).

Moreover, if constant solutions are not optimal, then there exists at least an inactivated arc $u(t) \equiv 1$. In particular, there are at most $2n$ switches in each regular domain.

In Fig. 2, we illustrate an example of evolution of the functions $\phi^{\pm}$ as functions of time along a non-singular normal trajectory. The latter shows how, due to the continuity of such function in a regular domain, the presence of an inactivated arc $u \equiv 1$ is required between two bang arcs.
## 2.5. Abnormal trajectories

Referring to Section 2.2.1, one can easily see that for abnormal trajectories admitting no singular arcs, the optimal control only takes the values $u \equiv u_{\min}$ or $u \equiv u_{\max}$ depending on the sign of $\phi = \langle F_2(x), p \rangle$ along the trajectory. It follows that there are at most $n-1$ switches between $u_{\min}$ and $u_{\max}$ in each regular domain, excluding the possibility of having inactivated arcs $u \equiv 1$, as illustrated on Fig. 2. As we shall see in the next section, abnormal extremals can in this case be seen as limit cases of normal extremals when the length of inactivated arcs tends to 0. This occurs because the endpoint at time $t_f > 0$ lies on the boundary of the accessibility set from the initial point $x_0$ at time $t_f$, which is a well-known geometric property of abnormal extremals (see, e.g., [30]). While, for the general case, the necessary conditions obtained through HMP do not give enough information about the admissibility of singular arcs along an abnormal extremal, they can be ruled out for the examples studied in the following section.

## 3. Applications

Previous works proved that simple two-dimensional systems are often capable of accounting for the main dynamical features of much larger GRNs [27,31]. The latter is due to the fact that most of the complex regulatory mechanisms observed in nature can be described by two main asymptotic patterns: multistability and oscillations. With this in mind, in this section, we focus on the two smallest GRNs motifs exhibiting bistability and damped oscillations, that can be both described by two-dimensional PWL models. While the biological objective is different in each case, in the control framework both problems can be reduced to inducing efficient state transitions. For certain values of the design parameter $\lambda$, the application of the HMP turns out to be informative enough to obtain explicit results concerning the optimal synthesis (i.e. a feedback control law). Both mathematical models obey the non-controlled two-dimensional dynamics defined in Filippov sense,

$$ \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^{\pm}(x_1, \theta_1), \end{cases} $$

where 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, and $s^{\pm}$ determines whether the system is a positive or negative feedback loop. The transcriptional thresholds delimit 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*} $$

which are defined as open sets in accordance with the HMP approach. We will see that controls are in general written as a concatenation of arcs characterized by a constant control. Thus, throughout this section, we resort to the notation $u_1 - u_2 - \cdots - u_n$ to imply that there exists a sequence $(t_k)_k$ of times such that $u(t) \equiv u_k$ for a.e. $t \in [t_k, t_{k+1}]$, for every $k$.

For each example, numerical simulations of both the open-loop and the resulting closed-loop dynamics are provided. Since the models here studied are described by simple two-dimensional dynamical systems, the uniqueness of the trajectories can be verified manually. However, in high-dimensional models of more complex GRNs, the use of the Filippov's formalism can produce multiple trajectories starting from the same initial condition, and thus it might be pertinent to verify the well-posedness of the resulting closed-loop system through more dedicated numerical methods. A more detailed discussion of this issue can be found in [32], where the authors propose a computational framework for GNRs based on reformulating PWL models as mixed complementarity systems. Additionally, since the controllers derived in this section are expressed as feedback laws, optimal trajectories can be obtained by simply simulating the closed-loop model (i.e., by injecting the feedback control $u(x)$ into the original open-loop model).

The theoretical results have also been verified using direct numerical methods through the Julia control-toolbox.<fnref n="2" /> This approach consists of discretizing the time variable to approximate the original OCP with a finite-dimensional optimization problem, that is subsequently solved with IPOPT (an interior-point optimizer). For that, the step functions $s^+$ and $s^-$ are regularized so that the discontinuous dynamics becomes continuous and nonlinear, contrary to the trajectories shown in the figures of this paper, which were obtained by a straight simulation of the PWL models. An important advantage of this method is that the optimization is not limited to predefined B-admissible trajectories. Thus, it can be helpful in performing preliminary studies of the system to determine suitable sequences of domains, as well as in confirming the choice of the sequence at the end of the theoretical study. The application of the direct methods is not included in this work, but for the sake of repeatability of the results, it can be found as an example of the Julia control-toolbox, where two regularization methods are compared (with sigmoid and exponential functions).<fnref n="3" />

<footnote n="2">² https://control-toolbox.org/.</footnote>

<footnote n="3">³ https://agustinyabo.github.io/PWLdynamics.jl/stable/.</footnote>
### 3.1. State transitions in a bistable switch

Consider two variables $x_1$ and $x_2$ which represent two genes mutually inhibiting each other. Then, the individual dynamics is

$$ \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 (3) $$

The domain $K = [0, k_1/\gamma_1] \times [0, k_2/\gamma_2]$ is forward invariant by the dynamics of Eq. (3), so that we consider only solutions evolving in $K$. Each regular domain $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), $$

and system (9) has two locally asymptotically stable steady states

$$ \phi_{10} = \left( \frac{k_1}{\gamma_1}, 0 \right) \in \bar{B}_{10}, $$

$$ \phi_{01} = \left( 0, \frac{k_2}{\gamma_2} \right) \in \bar{B}_{01}, $$

and an unstable Filippov equilibrium point at $(\theta_1, \theta_2)$. Fig. 3 illustrates the dynamical behavior of this non-controlled system. Following [7], we assume that both synthesis rates can be externally regulated, and so the controlled system, defined in Filippov sense, becomes

$$ \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}$. 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 Eq. (S) restricted to $\bar{B}_{00}$. In the coordinates $(x_1, x_2) \in B_{00}$, for $u \ge 1$, the separatrix $(S_u)$ can be written as the curve of equation

$$ 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}}, $$

which corresponds to the set of points in $B_{00}$ that belongs to the trajectory with constant $u(t) \equiv u$ for all $t$ starting in $x_2 = 0$ and reaching the point $(\theta_1, \theta_2)$. Using the latter, we define the regions

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

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

such that the domain $K$ is divided into

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

where $(\overline{S_u})^+$ and $(\overline{S_u})^-$ are forward invariant sets when Eq. (S) is driven by any constant control $u(t) \equiv u$ for every $t \ge 0$. Fig. 4 illustrates the role of these curves on the dynamic behavior of the system.

In boolean logic, the bistable switch is characterized by two possible states: $x \in B_{01}$ or $x \in B_{10}$, which are the regions of the state space associated to the stable steady states of the open-loop PWL dynamical system. The objective in this section is to induce transitions between these two regular domains. In particular, and without loss of generality, we focus on a transition from region $B_{10}$ to $B_{01}$, as the inverse problem is essentially equivalent due to the symmetry of the system. For that, we fix the sequence of regular domains $B = (B_{10}, B_{00}, B_{01})$, and restrict the optimal control problem to B-admissible trajectories, in the sense of Definition 3. Thus, we consider the boundary conditions

$$ x(0) = x_0 \in B_{10}, \quad x_2(t_f) = x_2^f \in \left( \theta_2, \frac{k_2}{\gamma_2} \right), $$

where $x_1(t_f)$ is not fixed since the production of gene $x_1$ is not controlled in $B_{01}$. For a given $\lambda \in (0, 1]$, our aim is to minimize $J_\lambda(u)$, as defined in Section 2.1. For feasibility of the control task (see [12]), assume the following.

**Assumption 2.** 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}. $$
<figure><img src="image_8.png" /><figcaption>Fig. 3. Streamplot of the open-loop bistable system (3) with parameters $k_1 = k_2 = 1$, $\gamma_1 = 1.1$, $\gamma_2 = 1.7$, $\theta_1 = 0.6$ and $\theta_2 = 0.4$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 4. Separatrix with parameters $k_1 = k_2 = 1$, $\gamma_1 = 0.25$, $\gamma_2 = 0.3$, $\theta_1 = 4$ and $\theta_2 = 3$, and $u_{\max} = 1.7$. For a fixed control $u$, if the state is in $(S_u)^+$ (respectively, $(S_u)^-$), then the state converges to $\phi_{01}$ (respectively, $\phi_{10}$).</figcaption></figure>

We start by recalling the following result for minimal time state transitions.

**Proposition 5** (Theorem 2 [12]). The minimal time strategy (i.e. $\lambda = 0$) is a feedback control law given by

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

Note that, for initial conditions in $(S_{u_{\max}})$, the latter control law yields optimal control functions with structure $u_{\min} - u_{\max}$, as shown in the example of Fig. 5. Naturally, since the term $|u - 1|$ related to control usage is not present in the cost function, no inactivated arcs are optimal. Using the simple structure of the Hamiltonian, one obtains the following result, which states that there is no singular arc along both normal (which is a direct consequence of Theorem 2) and abnormal trajectories. Its proof is given in Appendix C.

**Proposition 6.** No optimal trajectory $(x(t), u(t))$ can admit a singular arc. In particular, $u(t) \in \{u_{\min}, 1, u_{\max}\}$ for a.e. $t \in [0, t_f]$.

In the following theorem, we describe the possible structures of the optimal control for the general case.

**Theorem 7.** Along normal extremals, the optimal control is a feedback law

$$u(x) = \begin{cases} u_{\min} & \text{if } x \in (S_1)^{-}, \\ u_{\max} & \text{if } x \in B_{01} \text{ and } x_2 \ge x_2^s, \\ 1 & \text{otherwise,} \end{cases} \quad (5)$$

with

$$x_2^s(\lambda) = \frac{k_2}{\gamma_2} \frac{2\lambda - 1}{\lambda} \quad (6)$$
<figure><img src="image_10.png" /><figcaption>Fig. 5. Minimal-time trajectory with parameters $k_1 = k_2 = 1$, $\gamma_1 = 1.4$, $\gamma_2 = 1.6$, $\theta_1 = 0.6$ and $\theta_2 = 0.4$. The streamplot corresponds to the closed-loop control system resulting from the optimal control law (4).</figcaption></figure>

For abnormal extremals, the optimal control corresponds to the feedback law (4).

**Proof.** We start by proving the result for normal extremals ($p^0 \neq 0$), and, in accordance with the notations used in the proof of Proposition 6, we can define $(x(t))_{t \in [0, t_f]}$ (respectively, the adjoint state $(p(t))_{t \in [0, t_f]}$) as the concatenation of $(x^k(t))_{t \in [t_{k-1}, t_k]}$ for $k \in \{1, 2, 3\}$ (respectively, $(p^k(t))_{t \in [t_{k-1}, t_k]}$ for $k \in \{1, 2, 3\}$). In $B_{01}$ and $B_{10}$, one has at most one switch. The latter can be seen by the fact that the switching functions in $B_{01}$ are

$$ \phi^{-}(x, p, p^{0}) = k_{2}p_{2} + \lambda p^{0}, \qquad (7) $$

$$ \phi^{+}(x, p, p^{0}) = k_{2}p_{2} - \lambda p^{0} $$

and in $B_{10}$,

$$ \phi^{-}(x, p, p^{0}) = k_{1}p_{1} + \lambda p^{0}, $$

$$ \phi^{+}(x, p, p^{0}) = k_{1}p_{1} - \lambda p^{0}, $$

which can vanish at most once. As the feasibility of the trajectories implies that the control strategies start by $u \equiv u_{\min}$ in $B_{10}$ and end by $u \equiv 1$ or $u \equiv u_{\max}$ in $B_{01}$, we have that $u \equiv u_{\min}$ in $B_{10}$. Note that, since $x_1^0 < \frac{k_1}{\gamma_1}$, we cannot have $1 - u_{\min}$, as this would imply that during the 1 arc, the $x_1$-component of the trajectory is non-decreasing. Then, the length of the $u_{\min}$ arc (hence the cost $J_{\lambda}(u)$) would have to be larger. Moreover, we obtain a control structure $u \equiv 1 - u_{\max}$, $u \equiv 1$ or $u \equiv u_{\max}$ in $B_{01}$. The order of the arcs $u \equiv 1 - u_{\max}$ in $B_{01}$ is guaranteed by the fact that if there is a $u_{\max}$ arc in $B_{01}$, then $p_2$ is positive in $B_{01}$, and so the function $\phi^{-}$ is increasing in $B_{01}$. Now, suppose there is a switch in $B_{01}$ between the two arcs $u \equiv 1$ and $u \equiv u_{\max}$ at time $t_s$. Then, we have $\phi^{-}(t_s) = 0$, which implies that $p_2(t_s) = \lambda/k_2$. Using the fact that $x_1(t_f)$ is free, we have that $p_1(t_f) = 0$, and so $p_1(t) = 0$ for a.e. $t$ such that $x(t) \in B_{01}$. The Hamiltonian is

$$ H(x(t), p(t), p^0, u(t)) = \dot{x}_1(t)p_1(t) + \dot{x}_2(t)p_2(t) + p^0\lambda|1 - u(t)| + p^0(1 - \lambda) = 0. $$

Replacing the latter conditions, as well as $u(t) \ge 1$ and $p^0 = -1$, in the Hamiltonian evaluated at time $t = t_s$ yields (6), where $x_2^s = x_2(t_s)$. Corollary 1 applied with $n = 2$ implies that there are three possibilities in $B_{00}$: $u \equiv u_{\min} - 1 - u_{\max} - 1$, $u \equiv u_{\min} - 1 - u_{\max}$, or $u \equiv u_{\min} - 1$. The structure can be further explored by using the transversality conditions: using the fact that $x_1(t_f)$ is free, we have that $p_1(t_f) = 0$, and so $p_1(t) = 0$ for a.e. $t$ such that $x(t) \in B_{10}$ and the same holds in $B_{00}$, due to the continuity of the $p_1$ across these regions imposed in (T D). Thus, the switching functions in $B_{00}$ also have the form (7). Since we have $u \equiv u_{\min}$ in $(S_1)^-$, we obtain from Eq. (AD) that $p_2(t) < 0$ for a.e. $t$ such that $x(t) \in B_{00}$ and so no $u_{\max}$ arcs are allowed in $B_{00}$. Given that a unique $u \equiv u_{\min}$ is not admissible as the state never reaches $x_2 = \theta_2$, then a switch to $u \equiv 1$ is necessary. Now, consider a trajectory $(x(t))_{t \in [t_s, t_f]}$, where $t_s \ge 0$ is such that $x(t_s) \in (S_1)$, with control $u(t) \equiv u_{\min}$ for $t \in [t_s, t^*]$, then $u(t) \equiv 1$ for $t > t^*$, where the switching time between arcs is produced at time $t = t^* \ge t_s$. Let $t^* \mapsto T(t^*)$ be the function associating to $t^* \ge t_s$ the time at which we have $x_2(T(t^*)) = \theta_2$ with the latter control strategy. By computing the cost of such trajectory, we have

$$ J_{\lambda}(u) = \lambda(1 - u_{\min})(t^* - t_s) + (1 - \lambda)(T(t^*) - t_s). $$

Using that the function $T$ is increasing w.r.t. $t^*$ (as $\dot{x}_2|_{u=u_{\min}} < \dot{x}_2|_{u=1}$), we prove that $J_{\lambda}(u)$ is increasing w.r.t. $t^*$, and so the optimal switching time is $t^* = t_s$, which occurs exactly when the state $(x(t))_{t \in [0, t_f]}$ reaches the separatrix $(S_1)$.

Along abnormal extremals, i.e. $p^0 = 0$ as described in Section 2.5, there is at most one switch between $u_{\min}$ and $u_{\max}$ in $B_{00}$, while there is no switch between $u_{\min}$ and $u_{\max}$ in $B_{01}$ and $B_{10}$ since we have $\phi(x, p, p^0) = k_2p_2$ in $B_{01}$ and $\phi(x, p, p^0) = k_1p_1$ in $B_{01}$. By a similar computation of the cost $J_{\lambda}(u)$ to the one previously made for normal trajectories, it is possible to prove that the optimal switch occurs when the state reaches $(S_{u_{\max}})$. The feasibility of control strategies then proves the claimed structure. $\square$
<figure><img src="image_11.png" /><figcaption>Fig. 6. Optimal feedback control laws with parameters $k_1 = k_2 = 1$, $\gamma_1 = 1.4$, $\gamma_2 = 1.6$, $\theta_1 = 0.6$ and $\theta_2 = 0.4$. The colored regions indicate the value of the control to apply, and the streamplots illustrate the resulting closed-loop vector fields. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 7. Optimal trajectories for four different choices of $\lambda$ representing all possible control structures. Parameters values are $k_1 = k_2 = 1$, $\gamma_1 = 1.4$, $\gamma_2 = 1.6$, $\theta_1 = 0.6$ and $\theta_2 = 0.4$. The control bounds are $u_{\min} = 0.5$ and $u_{\max} = 1.5$.</figcaption></figure>

Fig. 6 illustrates both the abnormal and the normal feedback control laws in the $x_1x_2$-plane. In particular, we can see that for the “minimum fuel” state transitions given by $\lambda = 1$, there is no switch to the $u_{\max}$ arc in $B_{01}$, since $x_2^s(0) = \max(x_2) = k_2/\gamma_2$. This extreme case can be described by:

**Corollary 2.** The optimal control for the minimum fuel problem (i.e., $\lambda = 1$) is a feedback law of the form

$$u(x) = \begin{cases} u_{\min} & \text{if } x \in (S_1)^-, \\ 1 & \text{if } x \in (S_1)^+. \end{cases} \qquad (8)$$

Fig. 7 shows the optimal trajectories associated to four different values of $\lambda$, illustrating the different behaviors that can be obtained depending on the trade-off control usage/transition time. We conclude the study of the bistable switch by a description of these feedback control strategies, and a numerical analysis of the cost function.

**Remark 3** (Parametric Bifurcations w.r.t. the Trade-Off Parameter $\lambda$). As represented in Fig. 8, we identify four different optimal feedback regimes depending on the value of $\lambda \in (0, 1]$:

A Abnormal strategy (4),

B Normal strategy (5), with $u \equiv u_{\max}$ in $x \in B_{01}$,

C Normal strategy (5), with $x_2^s \in (\theta_2, \frac{k_2}{\gamma_2})$,

D Normal strategy (8).
<figure><img src="image_13.png" /><figcaption>Fig. 8. Plot of the cost of control (in gray), the final time (in orange), and the optimal switching time $t_s$ (in blue) as functions of $\lambda$ corresponding to the different optimal control strategies A, B, C and D described in Remark 3. The switching time $t_s$ corresponds to the time $t$ at which $x_2(t) = x_2^c(\lambda)$. Parameters values are $k_1 = k_2 = 1$, $\gamma_1 = 1.4$, $\gamma_2 = 1.6$, $\theta_1 = 0.6$ and $\theta_2 = 0.4$. The control bounds are $u_{\min} = 0.5$ and $u_{\max} = 1.5$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

The associated subintervals of $\lambda$ are delimited by the parameters $\lambda_a, \lambda_b, \lambda_c$, also illustrated in Fig. 8. The last two parameters can be easily computed from (6), while $\lambda_a$ requires an explicit computation and comparison of the costs associated to $\lambda = 0$ and $\lambda = \lambda_a$.

As a final remark, we can observe that every control switch along optimal trajectories occurs over one of the two separatrices. While the latter becomes useful in expressing the regulatory action as a feedback control, it is clearly a weak spot in terms of robustness. Indeed, if the switch is performed $\epsilon$ time before reaching the separatrix, the state returns to the initial regular domain, and the switch is not achieved. This phenomenon reveals another trade-off between optimality and robustness, which becomes critical in experimental settings under the presence of measurement noise inducing state uncertainty. While a thorough study of the problem is pertinent, as a preliminary solution, the latter can be compensated by delaying the control switch $\epsilon$ time after the separatrix is reached.

## 3.2. Sustained behavior in a damped genetic oscillator

The simplest GRN involving oscillatory behaviors can be modeled through two variables $x_1$ and $x_2$ that have opposite effects on each other: $x_1$ catalyzes the production of $x_2$, that in turn inhibits the production of $x_1$. Following the literature on the control of negative genetic feedback loops [33,34], we suppose that the system can be externally controlled by a chemical inducer that targets only one of the genes. Thus, we obtain the controlled PWL dynamical system

$$ \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 + k_2 s^{+}(x_1, \theta_1). \end{cases} \qquad (9) $$

The asymptotic behavior of the latter system has been extensively studied in the literature. In particular, in the open-loop case, it is a well known fact that constant control inputs cannot yield periodic trajectories, as they produce damped oscillatory behaviors towards its unique equilibrium:

**Lemma 2.** Under constant control functions $u \equiv u^*$, the state converges to the Filippov equilibrium point $(\theta_1, \theta_2)$ when $t \to \infty$.

**Proof.** The proof of this result can be found in [35, Theorem 2], where it is shown that for any constant control $u \equiv u^* > 0$, any trajectory converges asymptotically to $(\theta_1, \theta_2)$. $\square$

An example of this asymptotic behavior is shown in Fig. 9(a). This situation triggers a natural question: how can we use the external control to induce periodic trajectories? In this context, our objective is to obtain a cost-effective regulation law capable of transforming the non-controlled damped oscillations into sustained ones, while minimizing the cost $J_{\lambda}(u)$. To that end, we say an oscillation is sustained if, on each cycle, it passes by a given point in $\mathbb{R}^2$ which we define as the cycle point $(x_1^c, \theta_2) \in B_{10}$. Since the dynamics of (9) is not controlled in $B_{01}$ nor in $B_{11}$, it suffices to minimize $J_{\lambda}(u)$ among $B = (B_{00}, B_{10})$-admissible trajectories. Thus, we will first focus on the more general problem of reaching the cycle point from any initial condition in $B_{00}$, by fixing

$$ x(0) = x_0 \in B_{00}, \quad x(t_f) = (x_1^c, \theta_2) $$

for a free final time $t_f$. Once the general problem is grasped, we can tackle the more particular case of staying within the cycle, given by the more restrictive initial condition

$$ x(0) = (g(x_1^c), \theta_2), $$
<figure><img src="image_14.png" /><figcaption>Fig. 9. Open-loop behavior of system (9) (i.e. $u(t) \equiv 1$) converging asymptotically to $(\theta_1, \theta_2)$. The parameters are $k_1 = 2$, $k_2 = 4$, $\gamma_1 = 0.25$, $\gamma_2 = 0.3$, $\theta_1 = 4$ and $\theta_2 = 3$.</figcaption></figure>

where $g$ is the function defined as

$$g(x) = x \left[ \frac{\theta_2}{\theta_2 + \frac{k_2}{\gamma_2} \left[ \left( \frac{x}{\theta_1} \right)^{\frac{\gamma_2}{\gamma_1}} - 1 \right]} \right]^{\frac{\gamma_1}{\gamma_2}},$$

for $x \in [\theta_1, k_1/\gamma_1]$. For a given $x$, the value $g(x)$ corresponds to the $x_1$-coordinate of the point in the set $\{x \in \mathbb{R}^2 | x_2 = \theta_2, x_1 \le \theta_1\}$ reached by the trajectory starting in the cycle point, as illustrated in Fig. 9(b). The computation of the function $g$ is described in Appendix A. Similarly to what has been obtained in Proposition 6 for the bistable switch case, one has the following result concerning the absence of singular arcs:

**Proposition 8.** An optimal trajectory $(x(t), u(t))$ cannot admit singular arcs.

The proof, which follows exactly the same lines as the proof of Proposition 6, is based on a direct application of Theorem 1, with the choice $B = \{B_{00}, B_{10}\}$ and the controlled vector fields defined for $x = (x_1, x_2) \in \mathbb{R}^2$, $u \in [u_{\min}, u_{\max}]$, by

$$f_1(x_1, x_2, u) = \begin{pmatrix} -\gamma_1 x_1 + u k_1 \\ -\gamma_2 x_2 \end{pmatrix}, \quad f_2(x_1, x_2, u) = \begin{pmatrix} -\gamma_1 x_1 + u k_1 \\ -\gamma_2 x_2 + k_2 \end{pmatrix}.$$

In this setting, for $k \in \{1, 2\}$, the adjoint system is

$$\dot{p}_1^k = \gamma_1 p_1^k, \qquad (AD_1)$$

$$\dot{p}_2^k = \gamma_2 p_2^k,$$

and along an extremal, the maximized Hamiltonian $H$ satisfies

$$H(x(t), p(t), p^0, u(t)) = \dot{x}_1(t)p_1(t) + \dot{x}_2(t)p_2(t) + p^0 \lambda |1 - u(t)| + p^0 (1 - \lambda) = 0, \qquad (10)$$

for a.e. $t \in [0, t_f]$.

The following result on the structure of the optimal control can be obtained by direct application of the general theoretical result.

**Theorem 9.** Along normal extremals, an optimal control in a single regular domain $B_{00}$ or $B_{10}$ can be

* A single arc, which can be either a bang ($u \equiv u_{\min}$ or $u \equiv u_{\max}$) or an inactivated ($u \equiv 1$) arc,

* A concatenation of two arcs: a 1 arc followed by a bang arc ($u \equiv u_{\min}$ or $u \equiv u_{\max}$).

For abnormal extremals, the optimal control can only be of the form $u_{\min} - u_{\max}$ and $u_{\max} - u_{\min}$.

**Proof.** The proof simply holds from the fact that the switching functions in both regular domains

$$\phi^{-}(x, p, p^{0}) = k_{1}p_{1} + \lambda p^{0},$$

$$\phi^{+}(x, p, p^{0}) = k_{1}p_{1} - \lambda p^{0},$$
<figure><img src="image_15.png" /><figcaption>Fig. 10. Separatrices and nonadmissibility set of (9) for the cycle point $(x_1^c, \theta_2)$, with parameters $k_1 = 2$, $k_2 = 4$, $\gamma_1 = 0.25$, $\gamma_2 = 0.3$, $\theta_1 = 4$ and $\theta_2 = 3$; and control bounds $u_{\min} = 0.6$ and $u_{\max} = 1.4$.</figcaption></figure>

are monotone increasing or decreasing from the adjoint Equation (AD1) (depending on the sign of $p_1$), and so they can vanish at most once. Thus, depending on the value of $p_1$, they can yield either a single bang arc or a switch from a 1 arc to a bang arc. In order to prove the last claim, assume that $p^0 = 0$. Then, using the absence of singular arcs claimed by Proposition 8 it is easy to see that no switch may occur in a single regular domain, hence the control is equal to $u_{\min}$ or $u_{\max}$ in $B_{00}$ and $B_{10}$. Let us analyze all possible cases of signs by looking at (10) with $p^0 = 0$: if $p_2(t) > 0$ for a.e. $t \in [0, t_f]$, then $p_1(t) < 0$ for a.e. $t$ such that $x(t) \in B_{10}$, and $p_1(t) > 0$ for a.e. $t$ such that $x(t) \in B_{00}$, and so the control structure is $u_{\max} - u_{\min}$. Analogously, if $p_2(t) < 0$ for a.e. $t \in [0, t_f]$, then $p_1(t) > 0$ for a.e. $t$ such that $x(t) \in B_{10}$, and $p_1(t) < 0$ for a.e. $t$ such that $x(t) \in B_{00}$, and so the control structure is $u_{\min} - u_{\max}$. $\square$

As already established in Section 2.5, in optimal control problems with terminal constraints, the abnormal extremals can be very useful in characterizing the boundaries of the admissibility set. For that, we introduce some notation that will be important in describing admissible solutions: define the separatrix $(M_{u_1,u_2}) \in B_{00} \cup B_{10}$ as the stable manifold of the point $(x_1^c, \theta_2)$ under the control function

$$u(x) = \begin{cases} u_1 & \text{if } x \in B_{00}, \\ u_2 & \text{if } x \in B_{10}. \end{cases}$$

Analogously to the bistable switch case, the curve can be expressed as $x_2 = \delta(x_1, u_1, u_2)$, where the function $\delta$ can be computed by calculating the trajectory of the system under the given controllers, thus defining the sets

$$(M_{u_1,u_2})^+ = \{(x_1, x_2) \in \mathbb{R}^2 \mid 0 < x_1 < \frac{k_1}{\gamma_1}, \delta(x_1, u_1, u_2) < x_2 < \theta_2\},$$

$$(M_{u_1,u_2})^- = \{(x_1, x_2) \in \mathbb{R}^2 \mid 0 < x_2 < \theta_2, \delta(x_1, u_1, u_2) > x_2, 0 < x_1 < \frac{k_1}{\gamma_1}\}.$$

Then, based on the abnormal control structures $u_{\max} - u_{\min}$ and $u_{\min} - u_{\max}$, we use the separatrices $(M_{u_{\max},u_{\min}})$ and $(M_{u_{\min},u_{\max}})$ to define the domain of nonadmissibility of solutions

$$\mathcal{N}_a = (M_{u_{\min},u_{\max}})^+ \cup (M_{u_{\max},u_{\min}})^- \cup \{(x_1, x_2) \in B_{10} \mid x_1^c < x_1 < k_1/\gamma_1\}.$$

Both the separatrices and the nonadmissibility set are shown in Fig. 10. It is intuitive to see that, along an optimal trajectory we have $x(t) \notin \mathcal{N}_a$ for a.e. $t \in [0, t_f]$, but the trajectory can slide through the boundary $\partial \mathcal{N}_a$. Moreover, as seen in Fig. 10, the separatrix $(M_{u_{\max},u_{\min}})$ only spawns in $B_{10}$ for the current choice of parameters and control bounds, which is due to the inverse-time trajectory with $u \equiv u_{\max}$ that starts in $(x_1^c, \theta_2)$ not reaching $B_{00}$. However, this is not necessarily true in the more general case, for e.g., under weaker control actions (i.e., $u_{\min}$ and $u_{\max}$ closer to 1), or cycle points closer to the Filippov equilibrium.

So far, by application of the HMP, it was possible to reduce the original problem from an arbitrary constrained control $u \in [u_{\min}, u_{\min}]$ to a piecewise constant function composed of at most 4 arcs. In this approach, we further investigate the control structure by restraining the values of $\lambda$:

**Assumption 3.** $\lambda \in (0, 0.5]$.

Then, we see that the latter assumption allows one to reduce the space of admissible controllers by ruling out certain control structures:

**Lemma 3.** Under Assumption 3, the following holds along normal extremals:

* If $p_2(t) \ge 0$ for a.e. $t \in [0, t_f]$, there cannot be 1 arcs in $B_{00}$,
• If $p_2(t) \le 0$ for a.e. $t \in [0, t_f]$, there cannot be 1 arcs in $B_{10}$.

**Proof.** For the first claim, suppose $p_2(t) \ge 0$ for a.e. $t \in [0, t_f]$. Then, using (10) and Proposition 3, we get $p_1(t) > 0$ for a.e. $t$ such that $x(t) \in B_{00}$. By way of contradiction, suppose that there is a 1 arc in $B_{00}$ on the subinterval $[0, t_1]$, and thus, in that interval, we have $0 < p_1(t) < \lambda/k_1$ by Proposition 3. Condition (10) evaluated at $t \in [0, t_1]$ gives

$$(-\gamma_1 x_1(t) + k_1 u(t)) \underbrace{p_1(t)}_{<\lambda/k_1} + \dot{x}_2(t) p_2(t) - \lambda |u(t) - 1| - (1 - \lambda) = 0,$$

which then, using the fact that $|u(t) - 1| = 0$ for a.e. $t \in [0, t_1]$, yields the inequality

$$(-\gamma_1 x_1(t) + k_1) \frac{\lambda}{k_1} + \lambda - 1 > 0.$$

When solving for $x_1(t)$, and imposing $x_1(t) > 0$, we obtain

$$0 < x_1(t) < \frac{k_1}{\gamma_1} \frac{2\lambda - 1}{\lambda}.$$

which yields $\lambda > 0.5$, contradicting Assumption 3. For the second claim, suppose $p_2(t) \le 0$ for a.e. $t \in [0, t_f]$. Then, using (10) we get $p_1(t) > 0$ for a.e. $t$ such that $x(t) \in B_{10}$. Now, suppose that there is a 1 arc in $B_{10}$ on the subinterval $[t_s, t_2]$, where $t_s$ corresponds to the time at which $x_1(t_s) = \theta_1$, and thus, in that interval, one has $0 < p_1(t) < \lambda/k_1$. By a similar procedure as done in the first claim (replacing (10) in $t \in [t_s, t_2]$ to get an inequality, together with $|u(t) - 1| = 1$ and imposing $x_1(t) > \theta_1$), we can obtain $1 - \lambda > -\frac{\gamma_1}{k_1}\theta_1 + \lambda$, which implies $\lambda > 0.5$, again contradicting Assumption 3. □

Then, we can characterize the optimal control by investigating the admissible control structures:

**Proposition 10.** Under Assumption 3, the possible optimal control structures are

<table><thead><tr><th>Case</th><th>In B<sub>00</sub></th><th>In B<sub>10</sub></th></tr></thead><tbody><tr><td><b>A</b></td><td>1 - u<sub>min</sub></td><td>u<sub>max</sub></td></tr><tr><td><b>B</b></td><td>1 - u<sub>max</sub></td><td>u<sub>max</sub></td></tr><tr><td><b>C</b></td><td>u<sub>max</sub></td><td>1 - u<sub>max</sub></td></tr><tr><td><b>D</b></td><td>u<sub>max</sub></td><td>1 - u<sub>min</sub></td></tr></tbody></table>

where the 1 arcs can be of length zero depending on the initial condition.

**Proof.** The latter can be proved by studying the multiple combinations of signs of $p_1$ and $p_2$ in each regular domain:

• If $p_1(t) \le 0$ for a.e. $t$ such that $x(t) \in B_{00}$, then, using (10), we obtain $p_2(t) < 0$ for a.e. $t \in [0, t_f]$, and consequently $p_1(t) > 0$ for a.e. $t$ such that $x(t) \in B_{10}$. Again, Lemma 3 implies that the 1 arc in $B_{10}$ is not optimal, and so the structure of the optimal control is $1 - u_{min}$ in $B_{00}$ and then $u_{max}$ in $B_{10}$. The particular case $p_1(t) = 0$ in $B_{00}$ yields the abnormal control $u_{min} - u_{max}$.

• If $p_1(t) > 0$ for a.e. $t$ such that $x(t) \in B_{00}$, then the control can have 1 and $u_{max}$ arcs in $B_{00}$. Then:

- If $p_2(t) \le 0$ for a.e. $t \in [0, t_f]$, then $p_1(t) < 0$ for a.e. $t$ such that $x(t) \in B_{10}$, as $p_1(t) \le 0$ is forbidden by (10). Thus, only 1 and $u_{max}$ arcs are admissible in $B_{10}$. Lemma 3 states that the 1 arc in $B_{10}$ is not optimal. Then, the structure of the optimal control is $1 - u_{max}$ in $B_{00}$ and then $u_{max}$ in $B_{10}$. The particular case $p_2(t) = 0$ for a.e. $t$ yields the control $u_{max} - u_{max}$.

- If $p_2(t) > 0$ for a.e. $t \in [0, t_f]$, Lemma 3 states that the 1 arc in $B_{00}$ is not optimal. Then, there are three choices:

* If $p_1(t) > 0$ for a.e. $t$ such that $x(t) \in B_{10}$, then the structure of the optimal control is $u_{max}$ in $B_{00}$ and $1 - u_{max}$ in $B_{10}$.

* If $p_1(t) < 0$ for a.e. $t$ such that $x(t) \in B_{10}$, then the structure of the optimal control is $u_{max}$ in $B_{00}$ and $1 - u_{min}$ in $B_{10}$.

* If $p_1(t) = 0$ for a.e. $t$ such that $x(t) \in B_{10}$, then the structure of the optimal control is $u_{max}$ in $B_{00}$ and 1 in $B_{10}$. □

**Proposition 10** discards optimal control structures containing 1 arcs in both regions, as well as optimal control functions with a unique $u_{min}$ arc in a regular domain. A priori, the latter is not enough to fully characterize the optimal control: for each initial condition, one should compare the cost of each of the four control structures to obtain the optimal one. However, it is easy to see that each control structure is associated to a different region of $B_{00}$, and thus, for each initial condition $x_0 \in B_{00}$, there is only one admissible control structure, i.e., that can drive the state to the final cycle point $(x_0^c, \theta_2)$. The latter is illustrated in Fig. 11, that shows the four regions associated to each control structure, delimited by their respective separatrices.

Using this feature, it is then possible to express the optimal control as a feedback law.

**Theorem 11.** Under Assumption 3, the optimal control is defined for $x \in (B_{00} \cup B_{10}) \setminus \mathcal{N}_a$ by the feedback law

$$u(x) = \begin{cases} u_{\min} & \text{if } x \in (M_{u_{\max}, u_{\min}}) \text{ or } x \in (M_{u_{\min}, u_{\max}}) \cap B_{10}, \\ u_{\max} & \text{if } x \in (M_{u_{\max}, u_{\max}}) \text{ or } x \in (M_{u_{\max}, u_{\max}})^{-} \cap B_{00}, \\ 1 & \text{otherwise} \end{cases} \quad (11)$$
<figure><img src="image_16.png" /><figcaption>Fig. 11. Regions of $B_{00}$ linking each possible control structure (as defined in Proposition 10) to a set of initial conditions of the optimal control problem, delimited by the separatrices obtained for the cycle point $(x_1^c, \theta_2)$. The parameters are $k_1 = 2$, $k_2 = 3$, $\gamma_1 = 0.25$, $\gamma_2 = 0.3$, $\theta_1 = 4$ and $\theta_2 = 3$; and control bounds $u_{\min} = 0.6$ and $u_{\max} = 1.4$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 12. Optimal feedback control law for system (9) with parameters $k_1 = 2$, $k_2 = 4$, $\gamma_1 = 0.25$, $\gamma_2 = 0.3$, $\theta_1 = 4$ and $\theta_2 = 3$; and control bounds $u_{\min} = 0.6$ and $u_{\max} = 1.4$.</figcaption></figure>

For abnormal extremals, the optimal control structures are $u_{\max} - u_{\min}$ and $u_{\min} - u_{\max}$, which occur if and only if $x_0$ lies exactly over the separatrices $(M_{u_{\min},u_{\max}})$ and $(M_{u_{\min},u_{\max}})$, respectively, and so no arc is required to reach the cycle point.

The optimal feedback control is plotted in Fig. 12(a), while Fig. 12(b) shows different examples of optimal trajectories where the initial conditions have been chosen to produce the four possible control structures described in Proposition 10. As expected, every optimal trajectory uses an inactive arc ($u \equiv 1$) to reach a separatrix, and then slides over it until the cycle point.

For the particular case of starting from a point $(g(x_1^c), \theta_2)$ within the cycle, it is quite intuitive to see that the optimal control structures are either A or B depending on the choice of parameters and control bounds. Fig. 13 illustrates an example of a full cycle, where we see that the solution slides over $(M_{u_{\min},u_{\max}})$, thus having a $1 - u_{\min} - u_{\max}$ structure.

To conclude the study, we recall that the optimal control is not defined for $x \in \mathcal{N}_a$ due to the lack of feasibility of the control problem in that region. For the sake of completeness of the approach, we can formulate an extended feedback control law for every value of $x \in K$, capable of asymptotically driving the trajectories in $\mathcal{N}_a$ to the admissibility set, and subsequently applying the optimal feedback control (11) to reach the cycle point. Thus, for $x \in K$, we propose the extended optimality-based feedback control

$$u(x) = \begin{cases} u_{\min} & \text{if } x \in (M_{u_{\min}})^{-} \text{ or } x \in (M_{u_{\min},u_{\max}})^{+} \cap B_{10}, \\ u_{\max} & \text{if } x \in (M_{u_{\max},u_{\max}})^{+} \cap B_{10} \text{ or } x \in (M_{u_{\max},u_{\max}})^{-} \cap B_{00}, \\ 1 & \text{otherwise.} \end{cases} \quad (12)$$

whose closed-loop behavior is illustrated in Fig. 14(a). By definition, the latter control law is optimal for trajectories starting in $x \in (B_{00} \cup B_{10}) \setminus \mathcal{N}_a$. Moreover, as shown in Fig. 14(b), the trajectories starting in $\mathcal{N}_a \setminus ((\theta_1, \theta_2))$ under the proposed feedback control law converge to the cycle point $(x_1^c, \theta_2)$, and thus, to the optimal periodic trajectory, after a finite number of cycles.
<figure><img src="image_18.png" /><figcaption>**Fig. 13.** Optimal trajectory and control associated to a full cycle of system (9), with control structure A. The parameters are $k_1 = 2$, $k_2 = 3$, $\gamma_1 = 0.2$, $\gamma_2 = 0.3$, $\theta_1 = 4$ and $\theta_2 = 3$; and control bounds $u_{\min} = 0.4$ and $u_{\max} = 1.6$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 14.** Closed-loop behavior of system (9) under the extended feedback control strategy (12). The parameters are $k_1 = 2$, $k_2 = 3$, $\gamma_1 = 0.25$, $\gamma_2 = 0.3$, $\theta_1 = 4$ and $\theta_2 = 3$; and the control bounds are $u_{\min} = 0.6$ and $u_{\max} = 1.4$.</figcaption></figure>

## 4. Conclusion

In this paper, we developed a mathematical framework for inducing optimal state transitions in dynamical models of gene regulatory networks. For that, we considered a combination of the total time and $L^1$ cost of trajectories. The latter accounts for the main limitations of current practices in the external regulation of gene expression through chemical inducibles (such as IPTG), arising particularly in the field of synthetic biology and metabolic engineering. By integrating hybrid systems techniques with non-smooth optimal control methodologies, we managed to obtain general qualitative theoretical properties of the optimal control and the trajectories, that reduce the set of admissible controllers to piecewise constant functions composed of a limited number of arcs. In the second part of the paper, we focused on two application cases of PWL GRNs accounting for the main dynamical behaviors observed in biological regulatory mechanisms: oscillations and bistability. Both systems represent highly relevant examples from systems biology: the genetic toggle switch, a biosynthetic flip-flop device that acts as a biological memory unit through its bistability; and the genetic oscillator, capable of representing numerous naturally-evolved rhythmic behaviors in nature. For both cases, a thorough study of the problem reveals that the optimal control laws associated to each objective can be expressed in state-feedback form, and result in trajectories that slide over certain separatrices specific to the dynamics of each system. While these results remain mostly at the theoretical level, they can potentially provide guidance in designing practical closed-loop control schemes for experimental settings. We believe that the general framework could be a starting point for higher-dimensional studies: in that sense, we expect the illustrated bifurcations between abnormality and normality to be a generic feature of such problems. Moreover, we believe that exploiting the understanding of simpler two-dimensional motifs can help formulating cost-effective control laws for producing transitions in more complex GRNs of coupled synthetic devices.
Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Computation of the function g

Computing the g function can be done by explicitly calculating the trajectory of the system. We fix initial conditions $x_1(0) = x_1^c > \theta_1$ and $x_2(0) = \theta_2$, and so the solution for $t \ge 0$ of Eq. (9) in the first region $B_{11}$ is

$$ 
\begin{aligned}
x_1(t) &= x_1^c e^{-\gamma_1 t}, \\
x_2(t) &= \left(\theta_2 - \frac{k_2}{\gamma_2}\right) e^{-\gamma_2 t} + \frac{k_2}{\gamma_2}.
\end{aligned}
 $$

By defining $t_b \ge 0$ as the time such that $x_1(t_b) = \theta_1$, we obtain

$$ 
\begin{aligned}
x_1(t_b) &= \theta_1 = x_1^c e^{-\gamma_1 t_b}, \\
x_2(t_b) &= x_2^a = \left(\theta_2 - \frac{k_2}{\gamma_2}\right) e^{-\gamma_2 t_b} + \frac{k_2}{\gamma_2} = \theta_2 e^{-\gamma_2 t_b} + \frac{k_2}{\gamma_2} (1 - e^{-\gamma_2 t_b}),
\end{aligned}
 $$

which means that $t_b = \frac{1}{\gamma_1} \ln \frac{x_1^c}{\theta_1} > 0$, and thus $e^{\gamma_2 t_b} = \left( \frac{x_1^c}{\theta_1} \right)^{\frac{\gamma_2}{\gamma_1}}$. Then, in the second region $B_{01}$, we have

$$ x_2(t) = x_2^a e^{-\gamma_2(t-t_b)} = \left( \theta_2 + \frac{k_2}{\gamma_2} \left[ \left( \frac{x_1^c}{\theta_1} \right)^{\frac{\gamma_2}{\gamma_1}} - 1 \right] \right) e^{-\gamma_2 t}, $$

for every $t \ge t_b$. Defining $t_c \ge 0$ as the time at which $x_2(t_c) = \theta_2$, we get

$$ t_c = \frac{1}{\gamma_2} \ln \frac{\theta_2 + \frac{k_2}{\gamma_2} \left[ \left( \frac{x_1^c}{\theta_1} \right)^{\frac{\gamma_2}{\gamma_1}} - 1 \right]}{\theta_2}, $$

which means that

$$ x_1(t_c) = g(x_1^c) = x_1^c e^{-\gamma_1 t_c} = x_1^c \left[ \frac{\theta_2}{\theta_2 + \frac{k_2}{\gamma_2} \left[ \left( \frac{x_1^c}{\theta_1} \right)^{\frac{\gamma_2}{\gamma_1}} - 1 \right]} \right]^{\frac{\gamma_1}{\gamma_2}}, $$

for any $x_1^c \in [\theta_1, k_1/\gamma_1]$.

Appendix B. Proof of Theorem 2

**Proof.** In order to simplify the notations, we use $(x(t), u(t))$ to denote a given pair $(x^k(t), u^k(t))$ defined in a regular domain, for $t \in \Delta_k$. Assume that $u(t) > 1$ (respectively, $u(t) < 1$) for $t \in [t_1, t_2] \subset [0, t_f]$. Applying Theorem 1, the maximization condition provides that $u \equiv u_{\max}$ (respectively, $u \equiv u_{\min}$) outside the set $\phi^{-}(x(t), u(t), p(t), p^0) \equiv 0$ ($\phi^{+}(x(t), u(t), p(t), p^0) \equiv 0$, respectively). Without loss of generality, assume that $\phi^{-}(x(t), u(t), p(t), p^0) \equiv 0$ for a.e. $t \in [t_1, t_2]$, i.e. $\langle F_2(x(t)), p(t) \rangle + \lambda p^0 = 0$ for every $t \in [t_1, t_2]$. By successive differentiations w.r.t. $t$ (see Section 2.2.1), as the adjoint equation in Eq. (2) writes $\dot{p} = \Gamma p$ and the vector field $F_2(x(t))$ is constant for $t \in \Delta_k$, we obtain that $\langle F_2(x(t)), \Gamma^m p(t) \rangle \equiv 0$ for a.e. $t \in [t_1, t_2]$, for every $m \ge 1$. Let $q = \#\Xi_k$, define $\tilde{p}$ as the restriction of $p$ to the components belonging to $\Xi_k$, and denote the corresponding $q \times q$ reordered positive diagonal matrix built from $\Gamma$ by $\tilde{\Gamma}$. Thus we have $\langle K_2^k, \tilde{\Gamma}^m \tilde{p}(t) \rangle = 0$, for every $m \in \{1, \dots, q\}$, and $\langle K_2, \tilde{p}(t) \rangle \pm \lambda p^0 = 0$. As the $(\gamma_j)_j$ are pairwise distinct, we have that the family of vectors $(\tilde{\Gamma}K_2, \dots, \tilde{\Gamma}^q K_2)$ spans $\mathbb{R}^q$, and we can deduce $\tilde{p}(t) \equiv 0$, for every $t \in [t_1, t_2]$. Eq. (E) then provides $\tilde{p}(t) \equiv 0$ for $t \in \Delta_k$, and the second equality provides $\lambda p^0 = 0$, and hence $p^0 = 0$, so that $(x(t), u(t))_{t \in [0, t_f]}$ is an abnormal extremal, which is a contradiction. We can deduce that $u(t) > 1$ along a normal extremal for $[t_1, t_2] \subset [0, t_f]$ implies $u(t) \equiv u_{\max}$ (respectively, $u(t) < 1$ for $[t_1, t_2] \subset [0, t_f]$ implies $u(t) \equiv u_{\min}$), and the result follows. $\square$

Appendix C. Proof of Proposition 6

**Proof.** In order to apply Theorem 1, given the choice $B = \{B_{10}, B_{00}, B_{01}\}$, we set $\nu = 3$, and we define the vector fields, for $x = (x_1, x_2) \in \mathbb{R}^2$, $u \in [u_{\min}, u_{\max}]$, by

$$ f_1(x_1, x_2, u) = \begin{pmatrix} -\gamma_1 x_1 + u k_1 \\ -\gamma_2 x_2 \end{pmatrix}, \quad f_2(x_1, x_2, u) = \begin{pmatrix} -\gamma_1 x_1 + u k_1 \\ -\gamma_2 x_2 + u k_2 \end{pmatrix}, \quad f_3(x_1, x_2, u) = \begin{pmatrix} -\gamma_1 x_1 \\ -\gamma_2 x_2 + u k_2 \end{pmatrix}. $$
The times where changes of regular domains occur for the dynamics are denoted by $t_0 = 0 < t_1 < t_2$, and the final time is $t_3 = t_f$. We introduce the following functions $(\phi_j)_{j \in \{1, \dots, 10\}}$, which will guarantee the B-admissibility of the trajectories $x(t)$, which are solutions of Eq. (S). We define

$$z = (t_0, (t_1, x^1(t_0), x^1(t_1)), (t_2, x^2(t_1), x^2(t_2)), (t_3, x^3(t_2), x^3(t_3))),$$

where $t_0 = 0$ and $t_3 = t_f$, the latter time $t_f \ge 0$ being not fixed. In order to guarantee the B-admissibility and the continuity of the trajectory $(x(t))_{t \in [0, t_f]}$ at $t = t_1, t_2$, we define the functions

$$\left\{
\begin{aligned}
\phi_1(z) &= t_0, \\
\phi_2(z) &= x_1^1(t_0) - x_1^0, \\
\phi_3(z) &= x_2^1(t_0) - x_2^0, \\
\phi_4(z) &= x_1^1(t_1) - \theta_1, \\
\phi_5(z) &= x_1^2(t_1) - \theta_1, \\
\phi_6(z) &= x_2^1(t_1) - x_2^2(t_1), \\
\phi_7(z) &= x_2^2(t_2) - \theta_2, \\
\phi_8(z) &= x_2^2(t_2) - \theta_2, \\
\phi_9(z) &= x_1^2(t_2) - x_1^3(t_2), \\
\phi_{10}(z) &= x_2^3(t_f) - x_2^1.
\end{aligned}
\right.
$$

We can then define, for $\alpha = (\alpha_1, \dots, \alpha_{10}) \in \mathbb{R}^{10}$, the Lagrangian

$$L(z) = p^0 \left( \lambda \int_{t_0}^{t_v} |u-1|dt + (1-\lambda)(t_v - t_0) \right) + \sum_{j=1}^{10} \alpha_j \phi_j(z).$$

Applying Theorem 1 in this setting, Corollary 1 implies that the existence of a singular arc in a regular domain among $B_{10}, B_{00}, B_{01}$ implies $p^0 = 0$, i.e. the arc lies on an abnormal trajectory $(x(t), u(t))_{t \in [0, t_f]}$. Moreover, we have the following transversality and discontinuity conditions at times $t_0 = 0, t_1, t_2$ and $t_f$:

$$\left\{
\begin{aligned}
\alpha_1 &= 0, \\
p_1^1(0) &= \alpha_2, \\
p_2^1(0) &= \alpha_3, \\
p_1^1(t_1) &= -\alpha_4, \\
p_1^2(t_1) &= \alpha_5, \\
p_2^1(t_1) &= p_2^2(t_1) = \alpha_6, \\
p_2^2(t_2) &= -\alpha_7, \\
p_2^1(t_2) &= p_3^1(t_2) = \alpha_9, \\
p_3^2(t_2) &= -\alpha_8, \\
p_3^1(t_f) &= 0, \\
p_3^2(t_f) &= \alpha_{10}.
\end{aligned}
\right.
\qquad (TD)$$

Assuming $p^0 = 0$, for $k \in \{1, 2, 3\}$, the Hamiltonian $\mathcal{H}^k$ defined in Theorem 1 can be written as $\mathcal{H}^k = \mathcal{H}_0 + u^k \mathcal{H}_1^k$, with $u^k \in [u_{\min}, u_{\max}]$ where, for every $x^k = (x_1^k, x_2^k) \in \mathbb{R}^2$ and $p^k = (p_1^k, p_2^k)$,

$$\mathcal{H}_0(x^k, p^k, p^0) = -\gamma_1 x_1^k p_1^k - \gamma_2 x_2^k p_2^k,$$

$$\mathcal{H}_1^k(x^k, p^k) = \xi_1^k k_1 p_1^k + \xi_2^k k_2 p_2^k,$$

with $\xi_1^1 = 1$, $\xi_1^2 = 0$, $\xi_1^2 = 1$, and $\xi_2^2 = 1$, $\xi_3^1 = 0$, and $\xi_2^3 = 1$. In this setting, the adjoint State Equation (E) writes

$$\left\{
\begin{aligned}
\dot{p}_1^k &= \gamma_1 p_1^k, \\
\dot{p}_2^k &= \gamma_2 p_2^k,
\end{aligned}
\right.
\qquad (AD)$$

which is independent of $k \in \{1, 2\}$. For $k \in \{1, 2, 3\}$, singular arcs occur when the variables $(x^k(t), p^k(t), p^0, u^k(t))$ are extremal and satisfy $\mathcal{H}_1^k(x^k(t), \lambda^k(t)) = 0$, for every $t \in [T_1, T_2]$, where $0 \le T_1 < T_2 \le t_f$. Along such trajectories, the vanishing condition of the $k$th Hamiltonian $\mathcal{H}^k$ becomes

$$-\gamma_1 x_1^k(t) p_1^k(t) - \gamma_2 x_2^k(t) p_2^k(t) = 0, \qquad (V)$$

for every $t \in [T_1, T_2]$. Assuming that $\mathcal{H}_1^k(x^k(t), p^k(t)) = 0$ for $t \in [t_k, t_{k+1}]$ for some given $k \in \{1, 2, 3\}$, one can easily check that the vanishing condition on the Hamiltonian (V) together with discontinuity conditions (TD) imply that $(x(t), u(t))_{t \in [0, t_f]}$ is singular for
every $t \in [0, t_f]$, so that we obtain $\alpha_j = 0$ for every $j \in \{1, \dots, 10\}$. Hence the non-triviality condition of **Theorem 1** is violated, and the result follows. ☐

## Data availability

No data was used for the research described in the article.

## References

[1] N. Le Novere, Quantitative and logic modelling of molecular and gene networks, Nature Rev. Genet. 16 (3) (2015) 146-158.

[2] J. Ang, S. Bagh, B.P. Ingalls, D.R. McMillen, Considerations for using integral feedback control to construct a perfectly adapting synthetic gene network, J. Theoret. Biol. 266 (4) (2010) 723-738.

[3] L. Glass, S.A. Kauffman, The logical analysis of continuous, non-linear biochemical control networks, J. Theoret. Biol. 39 (1) (1973) 103-129.

[4] H. De Jong, J.-L. Gouzé, C. Hernandez, M. Page, T. Sari, J. Geiselmann, Qualitative simulation of genetic regulatory networks using piecewise-linear models, Bull. Math. Biol. 66 (2) (2004) 301-340.

[5] H. De Jong, J. Geiselmann, G. Batt, C. Hernandez, M. Page, Qualitative simulation of the initiation of sporulation in bacillus subtilis, Bull. Math. Biol. 66 (2004) 261-299.

[6] E. Farcot, J.-L. Gouzé, A mathematical framework for the control of piecewise-affine models of gene networks, Autom. 44 (9) (2008) 2326-2332.

[7] M. Chaves, J.-L. Gouzé, Exact control of genetic networks in a qualitative framework: the bistable switch example, Autom. 47 (6) (2011) 1105-1112.

[8] N. Augier, M. Chaves, J.-L. Gouzé, Qualitative control strategies for synchronization of bistable gene regulatory networks, IEEE Trans. Autom. Control 68 (2) (2022) 673-688.

[9] N. Augier, M. Chaves, J.-L. Gouzé, Weak synchronization and convergence in coupled genetic regulatory networks: Applications to damped oscillators and multistable circuits, Internat. J. Robust Nonlinear Control 33 (9) (2023) 4867-4892.

[10] P. Hillenbrand, G. Fritz, U. Gerland, Biological signal processing with a genetic toggle switch, PloS One 8 (7) (2013) e68345.

[11] N. Augier, A.G. Yabo, Time-optimal control of piecewise affine bistable gene-regulatory networks: preliminary results, IFAC-PapersOnLine 54 (5) (2021) 205-210, 7th IFAC Conference on Analysis and Design of Hybrid Systems ADHS 2021.

[12] N. Augier, A.G. Yabo, Time-optimal control of piecewise affine bistable gene-regulatory networks, Internat. J. Robust Nonlinear Control 33 (9) (2023) 4967-4988, http://dx.doi.org/10.1002/rnc.6012, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/rnc.6012 URL https://onlinelibrary.wiley.com/doi/abs/10.1002/rnc.6012.

[13] J.-B. Lugagne, S. Sosa Carrillo, M. Kirch, A. Köhler, G. Batt, P. Hersen, Balancing a genetic toggle switch by real-time feedback control and periodic forcing, Nat. Commun. 8 (1) (2017) 1-8.

[14] B. Shannon, C.G. Zamora-Chimal, L. Postiglione, D. Salzano, C.S. Grierson, L. Marucci, N.J. Savery, M. Di Bernardo, In vivo feedback control of an antithetic molecular-titration motif in escherichia coli using microfluidics, ACS Synth. Biology 9 (10) (2020) 2617-2624.

[15] L. Briand, G. Marcion, A. Kriznik, J.-M. Heydel, Y. Artur, C. Garrido, R. Seigneuric, F. Neiers, A self-inducible heterologous protein expression system in Escherichia coli, Sci. Rep. 6 (1) (2016) 33037.

[16] R.S. Donovan, C.W. Robinson, B. Glick, Optimizing inducer and culture conditions for expression of foreign proteins under the control of the lac promoter, J. Ind. Microbiol. 16 (1996) 145-154.

[17] F.C. Chittaro, L. Poggiolini, Singular extremals in $L^1$-optimal control problems: sufficient optimality conditions*, ESAIM: COCV 26 (2020) 99, http://dx.doi.org/10.1051/cocv/2020023.

[18] Z. Chen, J.-B. Caillau, Y. Chitour, $L^1$-minimization for mechanical systems, SIAM J. Control Optim. 54 (3) (2016) 1245-1265, http://dx.doi.org/10.1137/15M1013274, arXiv:https://doi.org/10.1137/15M1013274.

[19] A. Agrachev, B. Kazandjian, Optimal control for linear systems with $L^1$-norm cost, 2024, arXiv preprint arXiv:2404.07913.

[20] B. Berret, C. Darlot, F. Jean, T. Pozzo, C. Papaxanthis, J.P. Gauthier, The inactivation principle: mathematical solutions minimizing the absolute work and biological implications for the planning of arm movements, PLoS Comput. Biol. 4 (10) (2008) e1000194.

[21] H. Sussmann, A maximum principle for hybrid optimal control problems, in: Proceedings of the 38th IEEE Conference on Decision and Control (Cat. No.99CH36304), 1, 1999, pp. 425-430 vol.1, http://dx.doi.org/10.1109/CDC.1999.832814.

[22] F. Clarke, Necessary conditions, in: Functional Analysis, Calculus of Variations and Optimal Control, Springer London, London, 2013, pp. 435-471.

[23] T. Bayen, A. Bouali, L. Bourdin, The hybrid maximum principle for optimal control problems with spatially heterogeneous dynamics is a consequence of a pontryagin maximum principle for-local solutions, SIAM J. Control Optim. 62 (4) (2024) 2412-2432.

[24] T. Bayen, A. Bouali, L. Bourdin, O. Cots, Loss control regions in optimal control problems, J. Differential Equations 405 (2024) 359-397.

[25] A. Dmitruk, A. Kaganovich, The hybrid maximum principle is a consequence of pontryagin maximum principle, Systems Control Lett. 57 (11) (2008) 964-970, http://dx.doi.org/10.1016/j.sysconle.2008.05.006, URL https://www.sciencedirect.com/science/article/pii/S016769110800100X.

[26] T.S. Gardner, C.R. Cantor, J.J. Collins, Construction of a genetic toggle switch in Escherichia coli, Nature 403 (6767) (2000) 339-342.

[27] R. Edwards, S. Kim, P. Van Den Driessche, Control design for sustained oscillation in a two-gene regulatory network, J. Math. Biol. 62 (2011) 453-478.

[28] D. Ropers, H. De Jong, M. Page, D. Schneider, J. Geiselmann, Qualitative simulation of the carbon starvation response in Escherichia coli, Biosystems 84 (2) (2006) 124-152.

[29] T. Bayen, A. Bouali, L. Bourdin, O. Cots, On the reduction of a spatially hybrid optimal control problem into a temporally hybrid optimal control problem, in: Ivan Kupka's Legacy: A Tour Through Controlled Dynamics, Vol. 12, AIMS Applied Math., 2024, pp. 247-268.

[30] E. Trélat, Asymptotics of accessibility sets along an abnormal trajectory, ESAIM Control Optim. Calc. Var. 6 (2001) 387-414.

[31] R. Guantes, J.F. Poyatos, Dynamical principles of two-component genetic oscillators, PLoS Comput. Biol. 2 (3) (2006) e30.

[32] V. Acary, H. De Jong, B. Brogliato, Numerical simulation of piecewise-linear models of gene regulatory networks using complementarity systems, Phys. D: Nonlinear Phenom. 269 (2014) 103-119.

[33] I. Belgacem, J.-L. Gouzé, R. Edwards, Control of negative feedback loops in genetic networks, in: 2020 59th IEEE Conference on Decision and Control, CDC, IEEE, 2020, pp. 5098-5105.

[34] L. Chambon, I. Belgacem, J.-L. Gouzé, Qualitative control of undesired oscillations in a genetic negative feedback loop with uncertain measurements, Autom. 112 (2020) 108642.

[35] E. Farcot, J.-L. Gouzé, Periodic solutions of piecewise affine gene network models with non uniform decay rates: the case of a negative feedback loop, Acta Biotheor. 57 (4) (2009) 429-455.