Koopman Operators for Generalized Persistence of Excitation Conditions for Nonlinear Systems

TL;DR

This paper develops new identifiability conditions for nonlinear biological models using Koopman operators, linking data spectral properties to model predictability, and demonstrates their application with a synthetic gene circuit.

Contribution

It introduces a novel framework for nonlinear system identifiability based on Koopman operators and spectral data analysis, guiding experimental design for biological modeling.

Findings

Identifiability depends on dataset spectral rank in Koopman space.

Rank degeneracy leads to overfitting and poor model predictions.

Proposes criteria for designing experiments to ensure model predictability.

Abstract

It is hard to identify nonlinear biological models strictly from data, with results that are often sensitive to experimental conditions. Automated experimental workflows and liquid handling enables unprecedented throughput, as well as the capacity to generate extremely large datasets. We seek to develop generalized identifiability conditions for informing the design of automated experiments to discover predictive nonlinear biological models. For linear systems, identifiability is characterized by persistence of excitation conditions. For nonlinear systems, no such persistence of excitation conditions exist. We use the input-Koopman operator method to model nonlinear systems and derive identifiability conditions for open-loop systems initialized from a single initial condition. We show that nonlinear identifiability is intrinsically tied to the rank of a given dataset's power spectral density, transformed through the lifted Koopman observable space. We illustrate these identifiability conditions with a simulated synthetic gene circuit model, the repressilator. We illustrate how rank degeneracy in datasets results in overfitted nonlinear models of the repressilator, resulting in poor predictive accuracy. Our findings provide novel experimental design criteria for discovery of globally predictive nonlinear models of biological phenomena.

Equations

$${\mathbf{x}}_{t+1}=\mathbf{f}({\mathbf{x}}_t)$$

$$\mathcal{K}\psi(\bm{x}_{t})=\psi\circ\bm{f}(\bm{x})=\psi\big{(}\bm{f}(\bm{x}_{t})\big{)}=\psi(\bm{x}_{t+1}).$$

$$\begin{array}{rl}{\mathbf{x}}_{t+1}& =\mathbf{f}({\mathbf{x}}_t,{\mathbf{u}}_t)\\ {\mathbf{y}}_t& =\mathbf{h}({\mathbf{x}}_t)\end{array}$$

$$\mathcal{K}\psi ({\mathbf{x}}_t,{\mathbf{u}}_t)=\psi (\mathbf{f}({\mathbf{x}}_t,{\mathbf{u}}_t),{\mathbf{u}}_{t+1})=\psi ({\mathbf{x}}_{t+1},{\mathbf{u}}_{t+1})$$

$${\mathbf{\psi }}_x({\mathbf{x}}_{t+1})$$

$$\mathbf{y}$$

$$G_K(\text{z})={\mathbf{W}}_{\mathbf{h}}(\text{z}I-{\mathbf{K}}_x{)}^{-1}{\mathbf{K}}_{\mathbf{u}}$$

$${\mathbf{x}}_{t+1}$$

$${\mathbf{y}}_t$$

$${\mathbf{x}}_{t+1}$$

$${\mathbf{y}}_t$$

$${\mathbf{x}}_{t+1}$$

$${\mathbf{y}}_t$$

$$\mathbf{\psi }(\mathbf{f}(\mathbf{x}))=\mathbf{K}\mathbf{\psi }({\mathbf{x}}_t)$$

$$\begin{array}{rl}\mathbf{\psi }({\mathbf{x}}_{t+1})& =\mathbf{K}\mathbf{\psi }({\mathbf{x}}_t)+(\mathbf{\psi }({\mathbf{x}}_{t_0})-\mathbf{K}\mathbf{\psi }({\mathbf{x}}_t)){\delta }_{t,t_0-1}\end{array}$$

$$\begin{array}{rl}\mathbf{\psi }({\mathbf{x}}_{t+1})& =\mathbf{K}\mathbf{\psi }({\mathbf{x}}_t)+\mathbf{\phi }({\mathbf{u}}_t).\end{array}$$

$${\overline{\mathbf{R}}}_{\phi }(N):=\left[\begin{array}{ccc}{\mathbf{R}}_{\phi }(0)& \cdots & {\mathbf{R}}_{\phi }(N-1)\\ ⋮& \ddots & ⋮\\ {\mathbf{R}}_{\phi }(-(N-1))& \cdots & {\mathbf{R}}_{\phi }(0)\end{array}\right]$$

$$\bm{R}_{\varphi}(k)\coloneqq\mathbb{E}\big{[}\bm{\varphi}(\bm{u}_{t})\bm{\varphi}(\bm{u}_{t+k})^{T}\big{]},$$

$$\overline{\bm{R}}_{\varphi}(N)={\mathbb{E}}\bigg{\{}\begin{bmatrix}\bm{\varphi}(\bm{u}_{t+1})\\ \vdots\\ \bm{\varphi}(\bm{u}_{t+N})\end{bmatrix}\begin{bmatrix}\bm{\varphi}(\bm{u}_{t+1})^{T},\cdots,\bm{\varphi}(\bm{u}_{t+N})^{T}\end{bmatrix}\bigg{\}}$$

$$\begin{split}\overline{\bm{R}}_{\varphi}(N)&={\mathbb{E}}\big{[}\bm{v}\bm{v}^{T}\big{]}.\\ \end{split}$$

$$\begin{split}\bm{q}^{T}\overline{\bm{R}}_{\varphi}(N)\bm{q}&={\mathbb{E}}\big{[}\|\bm{q}^{T}\bm{v}\|_{2}^{2}\big{]}\\ &=\sum_{l,m=1}^{N}\bm{q}_{l}^{T}\bm{R}_{\varphi}(m-l)\bm{q}_{m}\\ &\geq 0\\ \end{split}$$

$$R_g(k)$$

$$\begin{gathered} \psi ({\mathbf{x}}_{t+1})=\left[\begin{array}{c}{\mathbf{x}}_{t+1}\\ \phi ({\mathbf{x}}_{t+1})\end{array}\right]=\left[\begin{array}{cc}{K}_{xx}& K_{x\phi }\\ K_{\phi x}& K_{\phi \phi }\end{array}\right]\left[\begin{array}{c}{\mathbf{x}}_t\\ \phi ({\mathbf{x}}_t)\end{array}\right]+ \\ {\delta }_{t,-1}(\mathbf{\psi }({\mathbf{x}}_0)-K\mathbf{\psi }(x_t)). \end{gathered}$$

$${\mathbf{S}}_{\phi }(\omega )=\sum _{k=-\infty }^{\infty }{\mathbf{R}}_{\phi }(k)e^{-ik\omega }$$

$$\begin{gathered} x_{t+1}=f({\mathbf{x}}_t)=K_{xx}{x}_t+K_{x\phi }\phi (x_t) \\ +{\delta }_{t,-1}({\mathbf{x}}_0-K_{xx}{\mathbf{x}}_t-K_{x\phi }\phi (x_t)) \end{gathered}$$

$${\mathbf{R}}_{\phi }(k)={\int }_{-\pi }^{\pi }{e}^{ik\omega }{\mathbf{S}}_{\phi }(\omega )d\omega .$$

$$\begin{split}\sum_{l,m=1}^{N}\bm{q}_{l}^{T}\big{[}\int_{-\pi}^{\pi}e^{i(m-l)\omega}\bm{S}_{\varphi}(\omega)d\omega\big{]}\bm{q}_{m}&\geq 0.\\ \end{split}$$

$$\begin{split}\int_{-\pi}^{\pi}\bigg{(}\sum_{l,m=1}^{N}\bm{q}_{l}^{T}e^{-il\omega}\bm{S}_{\varphi}(\omega)e^{im\omega}\bm{q}_{m}\bigg{)}d\omega&\geq 0,\\ \end{split}$$

$$\begin{split}\int_{-\pi}^{\pi}\bigg{(}\big{(}\sum_{l=1}^{N}e^{-il\omega}\bm{q}_{l}^{T}\big{)}\bm{S}_{\varphi}(\omega)\big{(}\sum_{m=1}^{N}e^{im\omega}\bm{q}_{m}\big{)}\bigg{)}d\omega&\geq 0.\end{split}$$

$$\int_{-\pi}^{\pi}\bigg{(}\bm{Q}^{*}(e^{i\omega})\bm{S}_{\varphi}(\omega)\bm{Q}(e^{i\omega})\bigg{)}d\omega\geq 0.$$

$$\sum _{k=1}^N{e}^{ik\omega }{\mathbf{q}}_k$$

$${\dot{m}}_{(i)}$$

Full text

I INTRODUCTION

Many physical systems exhibit phenomena with unknown governing dynamics. These systems can be high dimensional and partially modeled or completely unstudied. Self-assembling complex systems, biological networks, Internet-of-Things infrastructure models, smart cities, and social networks are all examples of dynamically evolving systems that frequently are represented by data. Identifying globally predictive models from data requires data collection strategies that capture all the modes of a system. For example, in biological network modeling and discovery, designing an informative set of experimental conditions can produce global biological models that capture multiple modes of dynamics, including invariant subspaces and multiple equilibria. Unfortunately, there are few identifiability metrics for quantifying the richness, or the informativity, of datasets of nonlinear systems. Further, a dataset may only elicit linear modes of a nonlinear system, even though nonlinear modes may lay dormant. In these scenarios, the accuracy of a model discovery algorithm is often confounded with the informativity or richness of a dataset used to train the model.

Nonlinear systems lack generalized criteria for quantifying the information content of a given dataset. Even if a nonlinear model is deemed globally identifiable, this may be under the assumption of continuous noise-free sampling (perfect data) of all states. In linear systems theory, the informativity of a dataset is characterized by its ability to persistently excite all the modes in the transfer function. These conditions, referred to as persistence of excitation conditions, prescribe rank requirements on either time-series or spectral signals. First introduced in [1], PE is “not (yet) consistently defined” [2]. There are two classical approaches to modeling persistence of excitation [1, 2] . Firstly, defining criteria or spectral properties of a control input such that it can elicit a response from the system, thereby uniquely specifying the frequency response model of the system ([3, 4]. Secondly, requiring that a control input that is non-zero in all channels, at least once throughout the time course ([5, 6]). Persistence of excitation conditions can also inform the appropriate construction of a universally exciting input signal, e.g. the construction of the Pseudo-Random Binary Sequence [7], which persistently excites all linear systems. It is thus desirable to establish a modeling framework to link classical results in linear system identification [3, 4, 6] for the treatment of nonlinear systems.

Recently, an emerging set of operator-theoretic tools have gained traction, centered on discovering linear representations of nonlinear dynamical systems in a lifted space of coordinates [8, 9, 10]. Originally derived for Hamiltonian systems [11], numerical [12] and theoretical [13] techniques for Koopman operator theory enable input [14, 15, 16, 17, 18, 19] and spectral modeling of nonlinear systems [9, 20], deep-learning based models of nonlinear phenomena [21, 22, 23], and study of chaotic and uncountable spectra arising in complex nonlinear phenomena [18, 24, 16, 8, 9]. In this work, we use the Koopman operator method (Sections II and III) to lift nonlinear systems into a linear space and thereby derive generalized persistence of excitation conditions (Section V), following a similar approach taken in linear systems theory. We illustrate how generalized PE conditions inform the design of initial conditions for simulated experiments on the three node repressilator (Section VI).

II KOOPMAN OPERATOR THEORY

In 1931, Koopman showed the existence of a coordinate transformation $\psi\in\mathbb{R}^{n_{L}}$ and a corresponding unitary operator ${\mathcal{K}}:\mathbb{R}^{n_{L}}\rightarrow\mathbb{R}^{n_{L}}$ for any Hamiltonian (non-dissipative) dynamical systems. He showed that the operator ${\cal K}$ and observable $\psi$ could be used to represent the time evolution of the underlying Hamiltonian system as a linear time-invariant system. This idea has been developed and generalized to other classes of nonlinear systems in recent years ([8, 9, 13], etc.), in the effort to find global representations instead of local approximations.

A discrete-time nonlinear dynamical system with state $\bm{x}_{t}\in\mathbb{R}^{n}$ at time $t\in\mathbb{N}$ under $\bm{f}:\mathcal{M}\rightarrow\mathcal{M}$ can be represented as:

[TABLE]

Then functions $\psi\triangleq\{\psi_{i}\}_{i=1}^{\infty}\in\mathcal{F},i\in\mathbb{Z}_{+}$ are called ”observables,” and represent a mapping from the state-space into a lifted set of coordinates. For example, $\psi$ may be a scalar observable comprised of nonlinearities, weighted combinations of the state, such as $\psi_{1}(\bm{x}_{t})=x_{1,t}\text{, }\psi_{2}(\bm{x}_{t})=x_{3,t}^{0.5},\psi_{3}(\bm{x}_{t})=x_{2,t}x_{n,t}^{2}$. The true Koopman observables are often approximated in numerical Koopman learning algorithms like extended DMD (eDMD), deepDMD, and Hankel dynamic mode decomposition using functions which define a generic basis on a Hilbert function space, e.g. radial basis functions (RBFs), Hermite polynomials, or a combination of such basis functions in deepDMD [21].

For an analytical function $f(\bm{x}_{t})$, we know from [14] that there exists a countably infinite or finite dimensional Koopman operator $\mathcal{K}$ that acts linearly on the observable $\psi:\mathcal{M}\rightarrow\mathbb{R}$ under function composition, satisfying the Koopman equation

[TABLE]

This equation states that the action of the Koopman operator on an observable function is equivalent to the action of the observable function on the vector field of the state.

From the above transient behaviour, one can see that this itself is a dynamical system in the space of observables whose time-evolution is governed by an infinite-dimensional operator that preserves the essence of (1). The transformation from (1) to (2) is frequently referred to as ”lifting” since the observable $\psi(\bm{x})\in\mathbb{R}^{n_{L}}$ usually has dimension $n_{L}\geq n.$ This is not always the case, but [23] shows that the observables define an expansion of the nonlinearities in the governing equations.

II-A Finite Dimensional Approximations

III Koopman Input Output Theory

Our study of the dynamical system (1), ultimately will require treating the initial condition as an input signal (via a Kronecker delta function) to the system. For clarity, we introduce the notion of an input-Koopman operator here, as well as the notion of the Koopman transfer function of a nonlinear system. Given a discrete-time nonlinear system with analytic vector field $\bm{f}(\bm{x},\bm{u})$ and control input $\bm{u}_{t}\in\mathbb{R}^{m}$, we write the dynamics of the system as:

[TABLE]

[16] showed that an observable $\psi:\mathcal{M}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{n_{L}}$ can lift the above system to $\mathcal{F}$ such that:

[TABLE]

For an exogenous memoryless input, [15] demonstrated that the above can be modified by splitting $\psi(\bm{x},\bm{u})$ into components $\psi_{x}(\bm{x})$ and $\psi_{u}(\bm{x},\bm{u})$, where $\psi_{x}(\bm{x})$ is a stacked vector valued observable of all scalar valued observable functions from $\psi(\bm{x},\bm{u})$ that do not depend on $\bm{u}$, and $\psi_{u}(\bm{x},\bm{u})$ is a stacked vector valued observable of all remaining terms. This results in the decomposed representation

[TABLE]

where $\tilde{\psi}_{u}({\bm{z}({\bm{u}_{t}})})\equiv\psi_{u}(\bm{x},\bm{u})$ and $\bm{z}_{t}$ is a stacked vector of all multivariate terms of $\bm{u}_{t}$ and $\bm{x}_{t}$. This is a representation of the nonlinear system dynamics that is linear in the lifted state observable $\psi_{x}(\bm{x})$ and the lifted input-state mixture observable $\psi_{u}(\bm{x},\bm{u})=\tilde{\psi_{u}}({\bm{z}_{t}})$. We thus can define the Koopman discrete time transfer function as

[TABLE]

where we have assumed that $\bm{y}_{t}=\psi_{x}({\bm{x}}_{t}).$ The zeros and the poles of the transfer function are defined in the usual manner, but with respect to the transformations on the state $\psi_{x}(\bm{x}_{t})$ and input $\tilde{\psi}_{u}(\bm{z}(\bm{u}_{t})))$. A continuous time analogue of the discrete-time Koopman operator is readily derived using the Koopman generator, rather than the Koopman operator. While the discrete-time formulation is chosen in this work, a similar treatment of the sequel can be used to develop persistence of excitation conditions for continuous-time nonlinear dynamical systems.

IV Problem Statement

In classical system identification theory, the design of an input or initial condition to guarantee identifiability is framed in terms of PE. One of the merits of a Koopman operator representation is that it provides a representation of a nonlinear system in linear coordinates. This allows us to rigorously formulate the question of how to design an initial condition that guarantees identifiability of the model. Specifically, we consider two problems:

Problem 1** (Identifiability of a Nonlinear System with Fixed Initial Conditions).**

Given a discrete time-invariant autonomous nonlinear system of the form

[TABLE]

determine if the model $\bm{f}(\bm{x})$ can be identified from the continuously sampled data stream $\bm{x}(t)$ and the set of initial conditions $X_{0}$.

This problem is difficult to solve, as it couples the problem of nonlinear function regression of $\bm{f}(\bm{x})$ and $\bm{h}(\bm{x})$ with the requirement of characterizing the richness of a set of fixed initial conditions. This problem is equivalent to the nonlinear state-space realization problem, given a single initial condition or set of initial conditions $X_{0}$.

The other variant of this problem is where the initial conditions are design parameters in an experiment and can be set by the user. This scenario is especially common in synthetic biology where experiments are conducted with varying initial conditions of concentrations, pH, etc. and hence is of interest from a design of experiment (DoE) standpoint.

Problem 2** (Identifiability of a Nonlinear System with Designed Initial Conditions).**

Given a continuous time-invariant autonomous nonlinear system of the form and the design choice of initial conditions $X_{0}$

[TABLE]

Find $X_{0}$ that guarantees the model $\bm{f}(\bm{x})$ can be identified from the continuously sampled data stream $\bm{x}(t)$.

This variant of the problem presumably has more degrees of freedom, namely $\dim(\bm{x}_{t_{0}})|X_{0})|=n|X_{0}|$ to be precise. However, the challenge is to relate identifiability of the model to the initial condition, which for an unknown $\bm{f}(\bm{x})$ is inherently difficult. We now use the input-Koopman framework derived in the previous sections to reformulate the problem with a linear Koopman representation. The recasting of the problem will permit extension of classical identifiability notions such as the design of PE or SR input signals or initial conditions.

V Persistence of Excitation Conditions for Nonlinear Systems with Koopman Operators

Our contribution in this paper is to formulate the problem of identifiability with fixed and designed initial conditions using the method of Koopman. Once we have formulated the problem in the Koopman operator theoretic framework, we derive computational certificates for PE for fixed or designed initial conditions. This leads to an algorithm for selection of the initial conditions of a dynamical system given $\bm{x}_{t_{0}}.$ In both scenarios, we treat the system’s initial condition as an input. Given a nonlinear discrete-time dynamic system of the form

[TABLE]

where $\bm{x}\in\mathbb{R}^{n}$ and $t\in\mathbb{Z},$ the Koopman equation for the corresponding system defines the action of an operator $\bm{K}$ on a vector valued observable $\bm{\psi}\in\mathbb{R}^{n_{L}}$ acting on the state $\bm{x}_{t}\in\mathbb{R}^{n}$, namely

[TABLE]

Then, (3) can be represented with initial state $\bm{\psi}(\bm{x}_{t_{0}})$ as an input:

[TABLE]

We thus model the initial condition as a Kronecker delta input to the dynamical system. Accordingly, we will abuse notation slightly and define the input of the system as $(\bm{\psi}(\bm{x}_{t_{0}})-\bm{K}\bm{\psi}(x_{t}))\delta_{t,t_{0}-1}\equiv\bm{\varphi}(\bm{u}_{t})$, which yields the classic input-state Koopman representation for a nonlinear system

[TABLE]

Notice that the state vector $\psi_{x}(\bm{x})\in\mathbb{R}^{n_{L}}$ is a vector observable function on which the Koopman operator acts as a linear operator to update the state. In essence, by defining an appropriate ’lifting’ or ’observable’ function, we can treat the problem as a classical linear identifiability problem, with additional caveats imposed by the presence of the nonlinearities in $\psi_{x}(\bm{x})$ and $\psi_{u}(\bm{x},\bm{u})$. As there are multiple definitions in the literature for PE [3, 6], we chose a definition that enables relating identifiability to the initial condition of a transfer function. We pose an extension of the definition from [3] which is that of SR from [6], but in the Koopman framework.

Definition 1** (Persistence of Excitation).**

A quasi-stationary discrete time observable $\bm{\varphi}(\bm{u}_{t})\in\mathbb{R}^{n_{L}}$ is said to be persistently exciting of order $N\in\mathbb{Z}_{+}$ if the covariance matrix $\overline{\bm{R}}_{\varphi}(N)$ is positive definite:

[TABLE]

where $\bm{R}_{\varphi}(k),k\in\mathbb{Z}$ is the auto-covariance of $\bm{\varphi}(\bm{u}_{t})$ formulated as:

[TABLE]

and ${\mathbb{E}}$ denotes the expectation operator. The above definition has been referred to as sufficient richness in [6] as well. Positive semi-definiteness of $\overline{\bm{R}}_{\varphi}(N)$ can easily be proved. By substituting the $R.H.S.$ of (10) in (9) we obtain:

[TABLE]

Defining \bm{v}\triangleq\big{[}\bm{\varphi}(\bm{u}_{t+1})^{T},\cdots,\bm{\varphi}(\bm{u}_{t+N})^{T}\big{]}^{T}, the above becomes:

[TABLE]

Let $\bm{q}=[\bm{q}_{1}^{T},\dots,\bm{q}_{k}^{T},\dots,\bm{q}_{N}^{T}]^{T},\bm{q}_{k}\in\mathbb{R}^{n_{L}}\setminus\{0\},$ we have

[TABLE]

Hence, $\overline{\bm{R}}_{\varphi}(N)$ is positive semi-definite. From the Herglotz Theorem, any function $R_{g}(k)$ defined on integers is positive semi-definite if and only if it has a Bochner representation (can be represented as the inverse Fourier transform) of a unique Spectral Measure $S_{g}(\nu)$ on a circle:

[TABLE]

If $R_{g}(k)$ is a scalar, $S_{g}(\nu)$ is its power spectrum as pointed out in [6] and [5]. Alternate definitions have been mentioned in literature but the following definition and result elucidates how PE of an initial condition for a dynamical system can be related to the Koopman transfer function.

Definition 2**.**

Given system (1), we say an initial condition $x_{0}$, treated as a Kronecker delta signal $\delta(\bm{x}_{0})$ is persistently exciting of Koopman-order $n_{L}$ if it is persistently exciting for a state-inclusive Koopman operator ${\cal K}\in\mathbb{R}^{n_{L}}$ of order $n_{L}$ satisfying

[TABLE]

Theorem 1**.**

The initial condition $x_{0}$ is persistently exciting for the nonlinear dynamical system (1) if and only if the Fourier transform of the auto-covariance matrix $\bm{R}_{\varphi}(k)$

[TABLE]

has $n_{L}$ distinct frequencies $\omega_{1},...,\omega_{n_{L}}$ where $\bm{S}_{\varphi}(\omega)$ does not vanish, i.e. $n_{L}$ positive spectral lines.

Proof.

Let the dimension of the lifted space and corresponding Koopman vector valued observable $\psi(\bm{x})$ be denoted as $n_{L}$. We write the initial condition as an input to the nonlinear dynamical system (1), of the form $\delta_{t}(\bm{x}_{0})$. We suppose that $\psi(\bm{x})$ is state inclusive [23]. This implies that:

[TABLE]

i.e. the PE of an input signal $\psi(\bm{x}_{0})$ for the system (13) implies the PE of the system (1), since the flow and vector field of system (13) is a projection of the flow and vector field of the system (1), respectively. It suffices to demonstrate the equivalence of the PE of $\psi(\bm{x}_{0})$, up to order $n_{L}$, to the linear independence of spectral lines for the power spectral density $\bm{S}_{\varphi}(\omega).$

We know that the spectral measure $\bm{S}_{\varphi}(\omega)$ of $\bm{R}_{\varphi}(k)$ would be a positive semi-definite, symmetric matrix of bounded measures, symmetric over all $\omega$ since the observable elements $\varphi(\bm{x})$ are all real valued. We have by the definition of the spectral power measure that

[TABLE]

Since the spectral measure is bounded almost everywhere, the monotone convergence theorem allows us to write from (12)

[TABLE]

After swapping the finite sum with limiting sums

[TABLE]

and distributing sums across the integrand, we obtain

[TABLE]

Defining filters $\bm{Q}(e^{i\omega})\triangleq\sum_{k=1}^{N}e^{ik\omega}\bm{q}_{k}$, the above becomes:

[TABLE]

Thus, positive definiteness of $\bm{S}_{\varphi}(\omega)$ only holds if and only if $\bm{R}_{\varphi}(k)$ is positive definite. But positive definiteness of $\bm{S}_{\varphi}(\omega)$ holds if and only if there are at least $n_{L}$ frequencies in which the integrand does not vanish, namely $n_{L}$ frequencies where

[TABLE]

is in the orthogonal complement of the null space of $\bm{S}_{\varphi}(\omega).$ This proves the result. ∎

From our simulation studies, we found that most initial conditions are PE up to order $n_{L}$, for their respective Koopman operator and dynamical system. If we consider the design problem of selecting an initial condition $\bm{x}_{0}$ such that $\psi(\bm{x}_{0})$ is PE, or a set of $\Psi(\bm{x}_{0})$ are PE of Koopman-order $n_{L}$, we can express the problem in terms of the positive definiteness of $\bm{R}_{\varphi}(N),$ adjusting the signal $\psi(\bm{x}_{0})$ or more generally the timing of the input to ensure the positive definiteness of the auto-covariance matrix. Alternatively, when working with a collection of initial condition signals $\psi(\bm{x}_{0})\in\Psi(\bm{x}_{0})$ it is straightforward to visualize the power spectrum using the transformed signal $\delta_{t}(\psi(\bm{x}_{0})\in\Psi(\bm{x}_{0})).$ Initial conditions can be selected or drawn randomly from the phase space until a suitable collection of initial conditions and $n_{L}$ spectral lines are identified.

VI Persistence of Excitation of Different Initial Conditions for a Repressilator Genetic Circuit Model

VI-A The Repressilator

The repressilator is a classical genetic circuit used in synthetic biology to implement circadian rhythms or synthetic oscillations. The architecture is that of a 3 node Goodwin oscillator, with three genes that produce proteins or mRNA that serve to repress the downstream or target gene’s function. Each gene represses its downstream target, with the final gene repressing the original gene to form a cycle of negative feedback. When the gain of the individual genes are balanced with respect to each other [25], the genetic circuit admits a limit cycle in the phase portrait and a single basin of attraction surrounding the origin.

There are many models of the repressilator, with varying degrees of complexity and intricacy to capture the underlying biophysical dynamics. We consider a simplified three dimensional model from the first experimental implementation of the repressilator [26], that captures the limit cycle and basin of attraction, to study the role of the initial condition in PE of the nonlinear system. Consider the model:

[TABLE]

An example set of commonly used initial concentrations is $1,0,0,0,0,0$ nM for LacI, $\lambda$-cI, TetR, mLacI, m$\lambda$-cI, and mTetR. We model the degradation and dilution rate of all proteins as a lump term with average kinetic rate $\delta=0.5$. Fig. 1(a) (and 2(a)) shows simulations (solid lines) of the repressilator from different initial conditions. The repressilator exhibits a strongly attracting limit cycle and a single unstable equilibrium point at the origin. Several initial conditions in the phase space mapped through the observable function have low gain, specifically those within $B_{1}(0)$ (unit ball in $\mathbb{R}^{6}$). These observable functions have low gain power spectra with spectral lines on the magnitude of numerical noise. We noted that these initial conditions, when mapped through higher order polynomials, lead to overfitting due to the vanishing of the signal in higher-order terms. To illustrate, $\bm{K}$ has been obtained using eDMD with third-order Hermite polynomials.

We considered training the repressilator model with initial conditions drawn from two different regions of the phase space. First, initial conditions within the unit disc centered at the unstable equilibrium point (solid lines in Fig. 1(b)) and secondly, initial conditions outside the unit disc centered at the same point (solid lines in Fig. 2(b)). We simulated using 6 initial conditions from one basin of attraction and evaluated test predictions using initial conditions from another basin. For example, we would train within the unit disc (Fig. 1(b)) or outside it (Fig. 2(b)) and evaluate prediction accuracy of the Koopman operator for trajectories initiated outside that basin of attraction (Fig. 1(c) and 2(c) respectively) using a norm based error between the predicted and simulated trajectories.

Notice the rank is greater for the power spectrum in Figure 1 than in Figure 2, which correlates with the failure to predict long-term global behavior in Figure 2. Interestingly, the rank of the power spectrum was not as high as the dimension of the list of dictionaries, indicating that the true Koopman observable space is of a lower dimension than the dimension of dictionary functions. Future work will investigate the iterative processes for identifying the minimal set of Koopman observables, their relationship to HankelDMD [27], as well as the formal design of automated biological experiments to ensure global predictive accuracy of discovered Koopman models.

VII ACKNOWLEDGMENTS

The authors would also like to thank Igor Mezic, Robert Egbert, Bassam Bamieh, Sai Pushpak, Sean Warnick, and Umesh Vaidya for stimulating conversations. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the Defense Advanced Research Projects Agency (DARPA), the Department of Defense, or the United States Government. This work was supported partially by a Defense Advanced Research Projects Agency (DARPA) Grant No. DEAC0576RL01830 and an Institute of Collaborative Biotechnologies Grant.