Stability Analysis of Parameter Varying Genetic Toggle Switches Using Koopman Operators

: The genetic toggle switch is a well known model in synthetic biology that represents the dynamic interactions between two genes that repress each other. The mathematical models for the genetic toggle switch that currently exist have been useful in describing circuit dynamics in rapidly dividing cells, assuming ﬁxed or time-invariant kinetic rates. There is a growing interest in being able to model and extend synthetic biological function for growth conditions such as stationary phase or during nutrient starvation. As cells transition from one growth phase to another, kinetic rates become time-varying parameters. In this paper, we propose a novel class of parameter varying nonlinear models that can be used to describe the dynamics of genetic circuits, including the toggle switch, as they transition from different phases of growth. We show that there exists unique solutions for this class of systems, as well as for a class of systems that incorporates the microbial phenomena of quorum sensing. Further, we show that the domain of these systems, which is the positive orthant, is positively invariant. We also showcase a theoretical control strategy for these systems that would grant asymptotic monostability of a desired ﬁxed point. We then take the general form of these systems and analyze their stability properties through the framework of time-varying Koopman operator theory. A necessary condition for asymptotic stability is also provided as well as a sufﬁcient condition for instability. A Koopman control strategy for the system is also proposed, as well as an analogous discrete time-varying Koopman framework for applications with regularly sampled measurements.


Introduction
Switching behavior in biological networks is of interest to many researchers across a variety of disciplines, such as epigenetics [1], synthetic biology [2], dynamical systems [3] and more [4][5][6]. The ability of complex biological networks to implement mixtures of switching behavior and nonlinear smooth system response make them a key area of focus in understanding how network structure, combined with system logic, are able to map structure to function [5][6][7][8]. The toggle switch is a special case of a two node Goodwin oscillator and thus provides the simplest biological realization of multi-stable behavior [9]. Finally, the toggle switch is present as a sub-network motif in many natural biological networks [4].
To date, most synthetic biological circuits are implemented in liquid suspension, assuming well stirred growth conditions and ample nutrient supply. Cells in this setting often divide exponentially, following a log-linear or log-phase growth curve. However, in every batch growth condition, bacterial cells will slowly deplete all nutrients in a growth medium, giving rise to late log-phase or early stationary phase dynamics [10]. In this setting, genetic circuits can be modeled as being subject to time-varying kinetic rates, as the cells are no longer in a state of dynamic, chemical equilibrium. Designing genetic circuits to function in nutrient-limited conditions is an open challenge.
Quorum sensing genes inside of living cells allow the cells to naturally sense when the cell population is becoming increasingly dense, saturating their growth environment, and potentially depleting resources. Quorum sensing genes activate in response to small homoserine lactone molecules, known as quorum signaling molecules, which are produced in cells containing these genes [11]. As the density of the cell population reaches saturation, the fraction of active quorum sensing genes begins to rise, which can act as a natural trigger for late log-phase or stationary phase biological processes. Several ordinary differential equation (ODE) models have been developed to model quorum sensing dynamics [12][13][14]. Nikolaev and Sontag modeled how a network of bacteria such as Pseudomonas fluorescens or E. coli use quorum sensing to synchronize gene expression according to toggle switches within each cell, thus showing how quorum sensing can synchronize switching behavior in a population of cells [15]. Here, we conceptually introduce the phenomenon of quorum sensing; in the sequel we will introduce a toggle switch model with time-varying kinetic rates that explicitly incorporates quorum sensing dynamics. Furthermore, we use quorum sensing in our models because we later implement a control strategy which requires the tracking of cell population density.
As a brief review, the classical genetic toggle switch [9] is modeled as follows: x = α 1 1 + y n y − δ 1 x (1) y = α 2 1 + x n x − δ 2 y where x ∈ R ≥0 and y ∈ R ≥0 represent the concentrations of repressor 1 (LacI) and repressor 2 (TetR); α 1 and α 2 denote the effective rate of synthesis of repressor 1 and 2; δ 1 > 0 and δ 2 > 0 are the decay rates of repressor 1 and 2 respectively; n x ≥ 1 and n y ≥ 1 represent the cooperativities of repression of promoter 1 (the Tet promoter) and 2 (the Lac promoter). This formulation can be represented using the standard mutual repression network model in Figure 1. A feature of the classical genetic toggle switch model is that the parameters are static, i.e., time-invariant. If we relax this assumption for the effective rates of synthesis and the decay rates, then we could have in our model that the kinetic parameters α 1 , α 2 , δ 1 , δ 2 now become time-varying functions, i.e., parameters that change over time. This generalization of the functional form of our system parameters, from time-invariant to time varying, provides a more accurate framework to model the genetic toggle switch in cells transitioning growth phases. Formally, the resulting mathematical model involving the time-varying parameters in question would then be a parameter varying nonlinear system (PVNS). Invoking this idea raises questions related to the stability of such a system as well as transient and steady state behavior. The analysis related to these PVNS toggle switch models is provided in this paper, as well as further findings pertaining to such a system with quorum sensing.
Koopman operator theory has been shown to provide insight into system behavior [16][17][18][19][20][21][22][23] through the analysis of the spectrum of a nonlinear system's Koopman operator (see Section 3 for a detailed introduction). Specifically, for biological applications, Ref. [24] showed that deep neural networks combined with Koopman operator methods could recover the governing equations of a glycolytic oscillator. Ref. [25] showed that Koopman operators could be used to formulate steady-state control of large-scale biological networks as a data-driven convex program, while Ref. [26] showed that cell growth could be modeled and predicted using Koopman operators that acted on temporal embeddings. Hasnain et al. [27] subsequently demonstrated that dynamic mode decomposition, a numerical method for estimating linear Koopman operators, could be used to quantify the impact of synthetic gene dynamics on a host organism. For observability and state estimation, Ref. [28] showed that Koopman operator methods could be used to determine sensor placement in sparsely measured gene networks and Ref. [26] showed that sensor fusion with Koopman operators enables more accurate identification of the underlying embedding space of Koopman observables.
Most recently, Johnson [29] showed that multi-stability of the genetic toggle switch can be successfully recapitulated with a special basis of Koopman observables that satisfied approximate closure. Nandanoori [30] showed that the spectrum of the genetic toggle switch's Koopman operator can be decomposed to mirror the invariant subspaces of the original toggle switch's phase space. In the literature described above, models of the toggle switch, or more generally, the biological network of interest assumed that cells were in a single growth phase, namely exponential growth phase. Motivated by the parametervarying nature of growth-phase transitions, here we define and extend the Koopman operator framework for the analysis of parameter varying nonlinear system models.

Parameter-Varying Toggle Switch (PVTS) System Model
The need for time-varying kinetic rates in the genetic toggle switch model comes from the current inability to model and extend synthetic biological function to growth conditions such as stationary phase or during nutrient starvation. With this quality in mind, the expression of the two genes of the toggle switch in the given organism can then be represented by the following parameter-varying nonlinear system.
where x ∈ R ≥0 , y ∈ R ≥0 represent the concentrations of repressor one and repressor two; α 1 (t) > 0 and α 2 (t) > 0 denote the effective rates of synthesis of repressor one and two. Note that we also require α 1 (t) and α 2 (t) to be bounded and to converge to positive constants; δ 1 (t) > 0 and δ 2 (t) > 0 are the self decay rates of repressor one, repressor two. Further, we note the requirement that δ 1 (t) and δ 2 (t) are bounded and converge to positive constants. The constants n x ≥ 1 and n y ≥ 1 represent the respective cooperativities of repression for promoter one and promoter 2; the Michaelis constants k x and k y are respective binding affinities. Similar to (1), we have that (x, y) ∈ R ≥0 × R ≥0 , i.e., the system is assumed to operate on the positive orthant. Given that α 1 (t), α 2 (t) are strictly positive functions over time, this is a parameter varying nonlinear system of the general form: where α(t) and δ are the vectors of time-varying parameters As a modelling choice, it makes sense to assume that our time-varying parameters are bounded because it would be physically improbable for there to be an unbounded rate of synthesis in α(t) or an unbounded rate of decay in δ(t). Furthermore, asserting these natural assumptions of boundedness will guarantee existence and uniqueness of solutions for this model (see Appendix A), which has already been established for the standard genetic toggle switch model.
This model is flexible in that α(t), δ(t) can be chosen to be sets of natural endogenous kinetic rates or to be sets of expression rate multipliers. As measurement technologies of kinetic rates in synthetic biology advance, it may be appropriate to select our time-varying parameters to be piecewise continuous functions, producing a piecewise smooth hybrid dynamical system, or to identify these parameter functions in a data-driven fashion. For the purposes of this paper, we will simply assume positivity, convergence to a positive constant, and boundedness.
Promoter leakiness has been modeled with the inclusion of time-invariant parameters which have been shown to affect stability properties of genetic toggle switches [31]. It is important to note that these parameters can be included in our models as multiplicative constants contained in the time-varying parameters and constant affine linear inputs. Using time-varying parameters to model resource usage is in contrast to other contemporary modeling approaches which use time-invariant parameters and a restructuring of the vector field [31].
Our models assume a negligible amount of cellular heterogeneity, as cellular heterogeneity would invoke the use of stochastic methods to model parameters [31]. This assumption holds when modeling bulk or average cell dynamics in batch or batch-fed reactor growth conditions; our framework thus is relevant for these physiological growth conditions. We remark that stochastic single-cell analysis with Koopman operator methods is beyond the scope of this paper and remains the subject of future studies. With a defined model class and physiologically motivated collection of assumptions, one important modeling criterion is that the generalized set of toggle switch equations should have unique solutions for the initial conditions on the domain. If we have an ODE of the form˙ x(t) = f ( x, t) where f ( x, t) is Lipschitz over some domain D in x with t ≥ t 0 , then we have existence and uniqueness of solutions to the initial value problem with initial conditions starting in D [32].
Since the classical genetic toggle switch has been shown to have unique solutions over its domain, then our more generalizing system should share this trait. This is given by Proposition 1 where the proof can be found in Appendix A.

Proposition 1.
Solutions to the PVTS with initial conditions x 0 starting in D = R ≥0 × R ≥0 exist and are unique for t ≥ 0.

Genetic Toggle Switch with Quorum Sensing
Cells with quorum sensing genes are able to regulate their gene expression based on the density of the cell population. We model the genetic toggle switch with quorum sensing (GTSQS) by incorporating the quorum sensing components of the genetic circuit using a 3rd state z that activates the first state x via a standard activator. The quorum sensing element u 3 will feed into this activator linearly, as will the respective inducers u 1 and u 2 of the two genes of the toggle switch x and y (see Figure 2 for a network diagram). The mathematical model of such a system is as follows: where x, y represent the concentrations of repressor 1 and repressor 2; z represents the quorum sensing activation response of the system that ultimately affects the concentrations x and y; α 1 (t) > 0 and α 2 (t) > 0 denote the effective rates of synthesis of repressor 1 and 2, respectively; α 3 (t) > 0 denotes the effective rate of the quorum sensing activator i.e., the rate at which quorum sensing affects gene expression; δ i (t) > 0, i ∈ {1, 2, 3} are the self decay rates of repressor 1, repressor 2, and the quorum sensing input, respectively; n x,y,z ≥ 1 represent the cooperativities of repression of promoter 1, promoter 2, and the quorum sensing input, respectively; k x,y,z > 0 are the respective binding affinities; (x, y, z) ∈ R ≥0 × R ≥0 × R ≥0 ; u 1 , u 2 are the inducers modeled as affine linear inputs; u 3 is the quorum sensing input. Note that we also require the time-varying parameters to be bounded functions of time which converge to positive constants.
In general, input signals to a biological network enter a nonlinear term, but the overall input signal may be generated by a Hill function [4]. Usually these input functions describe substrate saturation of a small molecule diffusing the cell membrane or a small molecule of mRNA designed to create interference, e.g., a CRISPRi gRNA designed to repress expression of a gene. In such a setting, there is always additive separability and so a generic input signal u(t) can be defined as the functional value of each of these nonlinear saturation functions. In short, this justifies the use of additive, linear input terms, even if the underlying substrate used to drive the control of the system may be subject to saturation. Further, for defined windows of input concentration, we can also argue that the inputs are approximately linear, when well within the windows of saturation.
This leads us to a GTSQS system of the form: where More generally, it is a parameter varying nonlinear system of the form:˙ where We present the system (7) below to showcase an example of GTSQS behavior.
For clarity, system (7) is the system (5) with the following parameters.
In this example, our time-varying parameters elicit gene expression responses we would expect in the context of nutrient starvation. On the time interval t ∈ [0, ∞), we see that our functions in α(t) go from 2 → 1 and our functions in δ(t) go from 1 3 → 1 2 as t goes from 0 → ∞. Intuitively, this means that the rates of synthesis are decreasing while the rates of decay are increasing. Furthermore, they satisfy our modeling requirements in that the functions are continuous, bounded, and exist in the limit as t → ∞.
From Figures 3 and 4 below, we can visualize the dynamics of example system (7), and get a brief intuition of the behavior of the GTSQS. We see in Figure 3 that the z = 0 plane is invariant, and there is one stable fixed point on this plane. We also see that for initial conditions with z > 0, go to one of two stable fixed points. Figure 4 shows that given a quorum sensing input u 3 = 0.01, we can take the initial conditions in the z = 0 plane to one of the other fixed points. In this way, our model captures how quorum sensing can affect log-phase and stationary phase gene expression.  . Simulated phase space of GTSQS with input u 1 = 0, u 2 = 0, u 3 = 0.01 in accordance with the example system in (7). Trajectories are plotted with initial conditions as blue empty circles and fixed points are seen as red filled circles.

Analysis of GTSQS Model
Here we provide mathematical analysis of our model over the domain D which is the first octant: The motivation for the choice of the domain of the GTSQS is that physically, it would be impossible for the state variables x, y, z (which represent quantities correlated to gene expression) to be negative in the sense of physical reality. Therefore, our model should not allow our state variables to be negative if we want the model to reflect biological phenomena. If one would like to experiment with possible negative states (outside the scope of this paper), then one should keep in mind that certain choices of n x , n y , n z may result in a vector field that have non-zero imaginary part, and inquiry into this topic would require tools from complex analysis.

Positively Invariant Domain
It can be trivially shown that for the genetic toggle switch (1) and for the PVTS (2) the first quadrant is positively invariant. Since this is the case, we require a similar result for our generalized model that is the GTSQS (4). We provide this result along with a proof that should serve as the blueprint for showing that the first quadrant is positively invariant for the systems (1) and (2). Proof. We want to show that the domain: We start our proof by assuming, for contradiction, that the system is not positively invariant, i.e., that for some (x 0 , y 0 , z 0 ) ∈ D, (x(t), y(t), z(t)) / ∈ D for some t ≥ 0. In order for this to be the case, the trajectory (x(t), y(t), z(t)) must have passed through the boundary of the first octant D B = {(0, y, z)} ∪ {(x, 0, z)} ∪ {(x, y, 0)} Plugging in these boundary coordinates to the dynamics of our system, it is clear that if x = 0, thenẋ ≥ 0; if y = 0, thenẏ ≥ 0; if z = 0, thenż ≥ 0 which implies that our trajectory could not have left D.

Remark 1.
We verify only for the u = 0 case because it is always mathematically possible to select some input u that would bring our state out of the positive orthant, even if such a control input would not be feasible in an actual genetic circuit.

Existence of Unique Solutions of GTSQS
In order to show that the GTSQS meets some basic criteria, and is in line with existing genetic toggle switch theory, we provide Proposition 2 to ensure that our model has unique solutions on our domain (proof found in Appendix B).

Proposition 2.
Solutions to the GTSQS with initial conditions x 0 starting in exist and are unique assuming bounded time-varying parameters where t ≥ 0.

Controllability of GTSQS
One of the properties to consider when designing a control system is controllability of the domain. In this section we see that the domain of the GTSQS is in fact controllable using the feedback linearization control strategy as follows: We can write the GTSQS as where whereū 1 ,ū 2 ,ū 3 are chosen to get desired steady state values. This implies that the system under the described control input has the desired GAS fixed point Since the desired fixed point can be chosen to be any point in our domain D, the system is controllable. Furthermore, since the controlled system is linear, the fixed point is unique, whereas the system with no input usually has more depending on the choice of parameters. This mathematical result from feedback linearization is not robust. We only showcase the rudimentary controller above to show that the system is in fact theoretically controllable. In order for this control law to work in practice, one must be completely accurate in capturing the nonlinearities and time-varying parameters in u(t), which is a non-trivial endeavor [32]. Hoping that this would work in a real genetic toggle switch is optimistic at best and a robust control strategy is required. The developments that have been made in the field of robust control allow us to synthesize theoretical controllers for this class of systems that give desired system behavior [33].

Stability Analysis of Parameter-Varying Nonlinear Systems with a Koopman Operator Approach
We now consider stability analysis of the genetic toggle switch, using a more general method that applies to a broader class of parameter-varying nonlinear systems. To begin, notice that the genetic toggle switch model with quorum sensing (2) can be written aṡ with the following assumptions : and f : R n → R n is a vector-valued analytic function, A(t) and ∆(t) are positive (all entries greater than zero for all t), bounded, continuous time-varying parameter matrices which converge to positive constant matrices in the limit as t → ∞.

Remark 2.
If one wishes to analyze the general form of the system in (9) or to look at a particular case of an example system by choosing parameters, take precaution in removing the assumptions established as you may not be able to guarantee that the initial value problem for (9) has unique solutions.
Parameter-varying nonlinear systems of this form have multiplicatively separable structure; the term A(t) f ( x) describes generalized birth processes in the creation of states while the term −∆(t) x describes generalized death processes. Both birth and death processes are attenuated or adjusted as a function of time, but the attenuation terms A(t) and ∆(t) may be independent and distinct for birth and death processes, respectively. As described previously, a strong motivation for systems with this structure is the presence of time-dependent kinetic rates in many biological reactions in living cells [10]. The assumptions previously stated allow us to arrive at conclusions that cater to these physical models of interest, while the most general form of results fall fairly directly from established theory on linear time-varying (LTV) systems [34].
Our first goal is to develop a generalized method to analyze the stability of this system. Direct analysis of a standard ansatz Lyapunov function such as V( x) = x T x can yield nice results, but in general, coming up with Lyapunov functions for stability is a nontrivial endeavor. We seek to extend the characterizations of stability from LTV theory to a broader class of nonlinear parameter varying systems that have the structure of Equation (9). To achieve this goal we will introduce and extend the method of Koopman to this class of parameter varying nonlinear systems because there has been much work done on the stability of nonlinear systems within this framework [19,[35][36][37][38]. First, we briefly introduce the Koopman operator representation of a time-invariant dynamical system.
Suppose that is a continuous time nonlinear dynamical system. There exists an infinite-dimensional linear operator U t that acts on a Hilbert function space G with invariant measure µ to satisfy the equation where K is known as a Koopman generator for the system (10), with corresponding Koopman observable g( x). In general, the Koopman generator admits an uncountable spectra, as it operates on an infinite dimensional function space. In special cases, say for when f ( x) is analytic, then the Koopman generator can be countably infinite dimensional or finite dimensional. The Koopman generator obtains its namesake in that it is used to generate a semigroup of flows on the state of the system, parameterized by the time variable t. In particular, when treating continuous time-invariant dynamical systems, we say that the Koopman generator generates a semigroup of Koopman operators [18].
Any element U t ∈ U t acts on the Koopman observable g( x 0 ), evaluated at some initial condition x 0 , to give g( x(t)) at time t ∈ [0, ∞). These elements U t are known as continuous-time Koopman operators that determine the flow of the dynamical system mapped through the Koopman observable g( x(t)), as opposed to the direct flow of the dynamical system in its natural state space [20].
The observable function g( x(t)) is the element of a function space G, which can be a function space of vector-valued, scalar-valued, matrix valued functions, etc. There are often multiple function spaces and multiple choices of Koopman observables that are consistent with an operator U t to satisfy the Koopman generator Equation (11), but a particular class of Koopman observables of interest are state-inclusive or state-recoverable Koopman observables [17,[39][40][41][42][43]. Such Koopman observables g( x(t)) satisfy the property that there exists some function P : G → R n such that x(t) = P( g( x(t))). An example of this is the class of Koopman observables where g( x) is a vector valued function of the form which is called a state inclusive Koopman observable. The nonlinear Koopman observables are held in ϕ( x) while the state of the system x is also used as an observable in g( x). The standard case of Koopman observables which is not state-inclusive is just g( x) = ϕ( x) [41]. The application problem often informs the appropriate function space for the Koopman observable.

State Inclusive Time-Invariant Koopman Learning
It is of interest to investigate the Koopman learning problem of Equation (9) which is a time-varying nonlinear system. Before we do this, let us cover the framework for state inclusive time-invariant Koopman Learning found in [39].
We begin with a nonlinear system of the forṁ and we choose our state-inclusive set of observables g( x) = x ϕ( x) such thaṫ and˙ which converts the system into the linear form: Note that here, we see that K xx is a n × n matrix K xϕ is a n × m matrix K ϕx is a m × n matrix K ϕϕ is a m × m matrix If m < ∞, then the matrix K is finite-dimensional for the chosen set of observables ϕ( x). In this case, we say that the basis of observables ϕ( x) satisfies the finite closure property. This case is a rarity in practice, and more than likely we choose a set of observables that gives an infinite-dimensional K matrix, but we proceed with an approximation of the system in the linear form by truncating the set of observables at some finite number m [39]. Lemma 1. If we have a nonlinear system that allows us to separate out the linear terms in the forṁ then there exists a time dependent Koopman Generator K(t) with state-inclusive observables Proof of Lemma 2. Choose our observables such that from (14), we havė and with (15), we arrive at the system in Equation (9) in the linear form

State Inclusive Time-Varying Koopman Learning
Previously, we have talked in great detail about the GTSQS, and now that we have presented it in the form of Equation (9), let us address state-inclusive Koopman Learning problem for this class of time-varying systems and obtain a time-dependent Koopman Generator [44].

Theorem 2.
If we have a nonlinear system that allows us to separate out the linear terms in the form of (9), then there exists a time dependent Koopman Generator K(t) with state-inclusive . (20) Proof. Starting with our system˙ is analytic, we can choose a set of observables, e.g., the Taylor polynomials, to give us˙ which allows us to arrive at the system in the linear time-varying (LTV) form for which we have our (n + m) × (n + m) time-varying Koopman generator K(t).

Spectral Properties of Parameter-Varying Koopman Generators
This LTV form˙ g( x) = K(t) g( x) allows us to extend the spectral results of LTV theory to Koopman models of the form in (20). From seeing that the system in Equation (9) gives a time-varying Koopman operator of the form (20), we begin with summarizing the "dynamic eigenvalue problem" [45] for this class of systems. We say that a dynamic eigenvalue λ(t) has an associated differentiable non-zero dynamic eigenvector u(t) if the eigenpair {λ(t), u(t)} satisfies in accordance with the classical LTV system theory [46]. The dynamic eigenvalue problem characterizes finding the dynamic eigenpairs {λ i (t), u i (t)} that satisfy Equation (21). This can be done classically with Wu's algorithm [46] stated as follows in the context of the time-varying matrix K(t) of size N := (n + m) 1.
Let the index variable j = 1 and let K j (t) = K(t) .

2.
Compute the eigenpairs of K j (t), denoted as {λ ij (t), u ij (t)}, i = 1, 2, ..., N and check if all of dynamic eigenvalues {λ ij (t)} are distinct. If they are distinct, then define the matrix S ij (t) which has the dynamic eigenvectors u ij (t) as its columns.
With this algorithm established as a way to solve the dynamic eigenvalue problem for our matrix K(t), let us also define the non-zero, differentiable vector (26) which is said to be the "mode vector" of K(t) associated with the eigenpair {λ i (t), u i (t)}.

Stability Properties of Parameter-Varying Koopman Generators
It is commonly known that one can use a Lyapunov function to guarantee the stability, asymptotic stability, exponential stability, or instability of a dynamical system [34]. One benefit of the Koopman form (20) is that we can use the previous section which summarizes the spectral properties of system (20) to write necessary and sufficient conditions for stability and asymptotic stability found in [46]. (ii) asymptotically stable if and only if that in addition to (i), we have that From Lemma 2, we see that we can deduce the stability of system (20) by solving the dynamic eigenvalue problem and checking to see if the conditions of Lemma 2 are satisfied for the system's mode vectors. This result applied to system (20) can be used to show the stability of the time-varying nonlinear system (9) without needing to find a Lyapunov function.
Since we have bounded time-varying coefficients for the LTV system (20), we can also use the necessary and sufficient condition for asymptotic stability through examining the "poles" of the system [47,48]. Let us briefly define the poles of the system, and summarize the condition for asymptotic stability. We have from (20) thaṫ Define a change of variables where L(t) is unitary, bounded, differentiable for t > t 0 , andL is bounded. L(t) is hereby a Lyapunov transformation, so it preserves the stability properties of the original system. L(t) brings the system (29) to˙ where (32) and the poles of the system are defined as the diagonal of the upper-triangular state matrix B(t). Let us denote this set of poles as P := {b ii (t)} ∀i ∈ {1, 2, ..., N}. Now we arrive at a lemma that allows us to determine the asymptotic stability of (29) from its poles.

Lemma 3. The system defined by (29) where K(t) is bounded is exponentially stable if and only if
Remark 3. Generally, coming up with a Lyapunov transformation L(t) that meets this requirements can be as difficult, but there exists algorithms that solve for this L(t); see [45].
Among these LTV results, there also exists methods of showing system stability that involve analysis on the system's isostables [37]. We previously mentioned the possibility of being able to use a Lyapunov function to show stability, which has been done in the Koopman framework by making use of the Koopman eigenfunctions [36,49,50]. In the context of our paper, the time-dependent eigenfunctions can be analyzed [44] and can also be used to formulate valid Lyapunov functions to show stability on a case-by-case basis.
When it comes to the design of a system, it would be desirable if the modelling portion of the design process was efficient. Since Wu's algorithm is computationally intensive, it would be desirable if we could easily check some conditions to ensure that the time-varying parameters chosen for the system have a chance at producing a stable system.
We know that the trace of K(t), denoted as tr(K(t)) is time-varying. Explicitly, we can write the trace as From Wu [46], we have the following proposition which provides a necessary condition for the asymptotic stability of a general LTV system.

Proposition 3. If the LTV system˙ x = A(t) x is asymptotically stable, then it is necessary for
We naturally extend this to the time-varying Koopman generator system with Lemma 4.
Now we use Lemma 4 to show that for Koopman systems of the form (20), we have a theorem that provides a simple way to determine if system (20) is not asymptotically stable. Theorem 3. If a PVNS with a Koopman representation of the form (20) with assumptions described as in (9) satisfies then the necessary condition of Lemma 4 is satisfied.
Proof. We apply Lemma 4 to our system (20) and find that if system (20) is asymptotically stable, then it is necessary for any t > t 0 that Recall the assumptions in (9) which state that ∆ ij (t) is bounded and ∆ ij (t) →∆ ij > 0 as t → ∞. So if we have that ∑∆ ij is some constant which is greater than the limit of tr(K ϕϕ (t)), then integral above goes to −∞.

Remark 4.
With this result, we can ensure the system (20) can be asymptotically stable. If the system fails to meet the condition of Theorem 3, then one can tune the parameters of ∆(t) to satisfy the condition. Note that just because a system in the form of (20) may satisfy the condition of Theorem 3, this does not mean the system is guaranteed to be stable because the condition in Theorem 3 only satisfies a necessary condition for asymptotic stability, not a sufficient one. Notice that the contrapositive of this theorem also provides an computational certificate to show the system is not asymptotically stable.
Note that it is possible to have asymptotic stability if but in this case, one must take the integral stated in Lemma 4 and ensure that it goes to −∞ to meet the necessary condition of Lemma 4.
We utilize a similar result from [46] in the proposition below to provide a sufficient condition for the instability of a system in the form of (20).
We now extend this to the time-varying Koopman generator system with Lemma 5.
then the system is unstable.

Proof.
We have the assumption that ∆ ij (t) is bounded and ∆ ij (t) →∆ ij > 0 as t → ∞. In addition, let us assume that ii < lim t→∞ tr(K ϕϕ (t)) .
Then we have that by our initial assumptions, we have that the above limit goes to ∞, which implies that the system (20) is unstable in accordance with Lemma 5.

Remark 5.
This theorem provides a sufficient condition for instability of a PVNS in the Koopman form of (20). Intuitively, one can understand the results from Theorem 3 and 4 as results on the steady state term ∆(t) which corresponds with the death processes in Equation (9). When the diagonal elements ∆ ii (t) are too small, the rates of decay or death in the biological network are too slow, which can lead to model instability.
The novelty of this section comes from the presented methodology of being able to take a PVNS with separable structure (9) to the time-varying Koopman form (20). In this form, the Koopman generator's diagonal steady state values determine whether the system can be asymptotically stable or not.

Nonlinear Feedback Control Design with Parameter-Varying Koopman Operators
It has been shown to be both possible and practical to implement control strategies within the realm of Koopman operator theory [49,[51][52][53]. We've shown in 2.3.5 that the GTSQS is mathematically controllable. Now, we wish to show that u(t) in the form of a negative observable feedback controller can stabilize systems of this form in the Koopman framework presented.
We take the time-varying system in Koopman form and linearly add in our controller u(t).
˙ ẋ We can rewrite the system by substituting giving us˙ Now in the same vein as in 2.3.5, we select u such that where we can select the blocks of the matrix M to give us stability. Depending on ∆(t) and A(t), it may be required that M xx or M xϕ vary with time in order to guarantee stability. Plugging this controller u(t) back into the original equations gives us the closed loop model which can be rewritten in the simplified forṁ where M(t) can be chosen to satisfy stability conditions.

Examples of Time-Varying Koopman Representations Being Used in Conjunction with Theorems 3 and 4 5.1. Example 1
Consider the system ẋ 1 This system is of the form (9) where So by letting ϕ( x) = x 2 1 we can write it in the form of (20) where Written more explicitly, we have: Now, if we attempt to design the system parameters such that we satisfy the condition of Theorem 3, then we must have that So assuming that we want our parameters to be bounded, the condition is satisfied as long as we choose µ(t) and λ(t) to converge to any positive constants.
From here, one can work through a method of one's choice (i.e., Wu's algorithm [46], Lyapunov function [54], poles via Lyapunov transform [47], etc.) to verify if a particular choice of µ(t) and λ(t) provide an asymptotically stable system.

Example 2
Consider the system where are all time-varying parameters are bounded and convergent. This system is of the form:˙ where So by letting we have brought the system into the Koopman form of (40) where Written more explicitly, we have that where Now, if we attempt to design the system parameters such that we satisfy the condition of Theorem 3, then we must have that So assuming that we want our parameters to be bounded, the condition is satisfied as long as we choose γ i (t) to converge to any positive constants. Also, note that the affine linear time-varying input u(t) permits us a limited fixed point selection. From here, one can work through a method of one's choice (i.e., Wu's algorithm [46], Lyapunov function [54], poles via Lyapunov transform [47], etc.) to verify if a particular choice of time-varying parameters provide an asymptotically stable system. Remark 6. The implication of satisfying the condition of Theorem 3 has been shown using the methodology described in this paper. If one wishes to satisfy the condition of Theorem 4 for the two examples presented, then flip the inequality symbol in (43) and (47) and choose appropriate parameters that cause the resulting systems to be unstable.

Discrete Time Analog
Since the Koopman operator framework also extends to discrete-time systems [18,55], we show a brief extension of the time-varying Koopman operator now in discrete-time. Suppose we have a PVNS in discrete-time of the following form: Then, similarly to the process in Section 4, by choosing the state-inclusive set of such that x(t n ) + A(t n )K xϕ ϕ( x(t n )) (49) and we obtain the discrete-time-varying system in the following form: Remark 7. By discretizing (9) with the Forward-Euler method, we can arrive at the discrete- where h n := (t n+1 − t n ). We can deduce the stability properties for the discrete time systems (50) and (51) by using the well established results on the stability of discrete time-varying linear systems found in [56]. In applications, the gene expression data from a genetic toggle switch comes as a sequence of vectors. The new form (51) explicitly shows how the current instance of the parameter functions affect the next snapshot of data observed. This is in contrast to purely data-driven approaches such as dynamic mode decomposition, where the effect of time-varying parameters and the evolution of the system states are not easy to deconvolve.

Conclusions
In this paper, we extended the classical genetic toggle switch model to a parameter varying toggle switch that maintained the key characteristics of the original model while allowing for the kinetic rates of synthesis to vary with time. Next, we provided the model of the genetic toggle switch with quorum sensing which includes the one sided activation of a gene due to quorum sensing, as well as affine linear inducer inputs. We then showed that the genetic toggle switch with quorum sensing is controllable using feedback linearization. In addition to the genetic model provided in the paper, we showed how we can use Koopman operator theory to take a parameter varying nonlinear system with a particular separable structure to the form of a time-varying Koopman generator equation. Further, we showed how to determine the stability of the system in its Koopman generator form using theorems derived from classical linear time-varying system theory. Within this Koopman framework, we exhibited a general control strategy for the time-varying Koopman system analogous to the one provided for the genetic toggle switch with quorum sensing. Finally, we utilized the time-varying Koopman framework in the discrete-time setting, showing that if we had a parameter varying nonlinear discrete-time system with a particular separable structure, then we can use Koopman operator theory to take this system to the form of a time-varying Koopman operator equation.  − α 2 (t)n x (x/k x ) (nx −1) We have the assumption that 0 < α 1 (t) ≤ᾱ 1 for some constantᾱ 1 , 0 < α 2 (t) ≤ᾱ 2 for some constantᾱ 2 , 0 < δ 1 (t) ≤ δ 1 for some constant δ 1 , 0 < δ 2 (t) ≤ δ 2 for some constant δ 2 .
By examining the off-diagonal entries of J, we see that these functions are bounded because the degree of the leading power term in the denominator is greater than that of the numerator in both terms. With this argument, we have that α 1 (t)n y (y/k y ) (n y −1) ((y/k y ) n y + 1) 2 ≤ᾱ 1 n y , − α 2 (t)n x (x/k x ) (n x −1) ((x/k x ) n x + 1) 2 ≤ᾱ 2 n x .

Appendix B
Here we prove Proposition 2.
Proof. We will show that the Jacobian of the dynamics of the GTSQS is bounded over the subset of the domain where z > 0. Later in the proof, we show that the vector field of the GTSQS is Lipschitz for z = 0, thus showing that the entire domain is Lipschitz. We can partition the domain in this way because the z = 0 plane is invariant under the vector field of the GTSQS. We will continue by bounding the infinity norm of the Jacobian. The Jacobian is: −α 1 (t)( z kz ) nz ( 1 ky ) ny n y y ny −1 (1+( y ky ) ny +( z kz ) nz ) 2 α 1 (t)n z z(1+( y ky ) ny +( z kz ) nz ) 2 −α 2 (t)n x x nx −1 ( 1 kx ) nx (1+( x kx ) nx ) 2 −δ 2 (t) 0 0 0 By ensuring the boundedness of each component of our Jacobian matrix, we verify the boundedness of the infinity norm of the Jacobian. The zero components, of course, are trivially bounded. We are assuming bounded time-varying parameters, so we explicitly write these bounds as 0 < α 1 (t) ≤ᾱ 1 for some constantᾱ 1 , 0 < α 2 (t) ≤ᾱ 2 for some constantᾱ 2 , 0 < α 3 (t) ≤ᾱ 3 for some constantᾱ 3 , 0 < δ 1 (t) ≤ δ 1 for some constant δ 1 , 0 < δ 2 (t) ≤ δ 2 for some constant δ 2 , 0 < δ 3 (t) ≤ δ 3 for some constant δ 3 .
Let us examine the following function: Next, consider the function: We have that the power terms in the denominator are strictly greater than that of the numerator, thus By similar argument, we have that |J 12 | ≤ᾱ 1 n y ( 1 k y ) n y .
This function is undefined for z = 0, so in order for this function to be bounded, we must require that z > 0. With this, we have |J 13 | ≤ᾱn z .
Now we show that the vector field is Lipschitz for z = 0 with Lipschitz constant L z=0 = max{ δ 1 +ᾱ 2 k x , δ 2 } .
From (4), we see that for z = 0, u = 0, we have that Let x, x ∈ D z=0 . Then we want to show that for some constant L z=0 , we have the following: Let us proceed with n x = 1 for simplicity and without loss of generality, and the norm is the 1-norm.