Prediction of fitness in bacteria with causal jump dynamic mode decomposition

One of the most fundamental processes in life is the ability to replicate and pass on hereditary material [1]. From viral particles to bacteria to mammalian cells, cell division is fundamental to growth, maintenance of physiological health, and intrinsically tied to the notion of senescence [2].

The mechanisms for controlling growth in organisms are determined by metabolic networks [3, 4], namely their topological structure and parametric realization. Known metabolic networks in well studied model organisms such as E. coli [5] and S. cerevisiae [6, 7] have given rise to predictive models that translate environmental activity to metabolic network state, and ultimately to predictions of growth rate. For canonical biological model systems, these models have been highly accurate in predicting growth rate and found utility in industrial microbiology applications, e.g. in the design of bioreactors or informing best practices in food safety.

For many biological life forms, relatively little is known about their metabolic network or structure. This is especially the case when developing bioengineering tools in novel host microbes [8, 9]. For new organisms, canonical metabolic networks are lacking and often obtained through a process of sequence alignment and comparative analysis with existing metabolic network models in relative species. However, many novel strains do not exhibit significant similarity, and even in the case of sequence similarity, small mutations can lead to dramatically different growth phenotypes, e.g. growth of non-pathogenic soil strains[10, 11] versus pathogenic counterparts [12]. The absence of predictive cross-species models, as well as the inability to predict growth phenotype wholly from sequence data, motivates the need for data-driven methods to accelerate the discovery of metabolic models and growth rate prediction models.

Due to advances in high-throughput experimental techniques, it is relatively easy to characterize growth rates as a function of exposure to environment. Liquid and acoustic-liquid handling robotics enables interrogation of thousands of growth conditions in a single microtiter plate, which in turn opens the door for using data-driven approaches [13] to predict growth rate as a function of environmental state. Is it possible to accurately predict the growth rate of a microbe, entirely from the chemical composition and environmental parameters of its growth condition? In this paper we explore a data-driven operator theoretic approach that utilizes microtiter plate reader data, and more generally multi-variate time-series data, to develop predictive models of growth rate in Pseudomonas putida, a broadly used strain for commercial bioreactors and a target workhorse for tractable genetic engineering [14].

A broadly successful class of data-driven modeling approaches stem from the study of Koopman operators, a mathematical construct for representing the time-evolution of nonlinear dynamical systems. In Koopman operator theory, the time-evolution of a nonlinear system is defined on a function space, acting on the original state of the system. In this function space, known the observables space, the Koopman operator is a linear operator, enabling spectral analysis, the decomposition of eigenspaces, and study of nonlinear structure [15]. The Koopman operator framework has been developed for continuous [16] and discrete time systems [17, 18], for open-loop [17] and input-controlled [19, 20] dynamic systems. Thus, Koopman operators present a powerful framework for analyzing the behavior of nonlinear systems, including predicting how experimental conditions regulate growth dynamics.

Many numerical methods for identifying Koopman operators from data have been developed in the last two decades [21, 22, 23, 24, 25, 26, 27, 28]. The most common approach is to use dynamic mode decomposition (DMD), which models nonlinear dynamics via an approximate local linear expansion [25]. In [17] an extended dictionary of basis functions with universal function approximation properties is used to discover an approximation of the lifting map or observables. These techniques suffer from combinatorial explosion, which generally has prohibited analysis of high-dimensional nonlinear systems [29]. The most recent developments in the field of DMD integrate established advances in deep learning with DMD [27, 30, 31, 32] where the deep neural networks have the capacity to approximate exponentially many distinct observable functions. Recent work has shown deep Koopman learning algorithms can be extended to synthesize controllers for systems subject to uncertainty [33], suggesting that deep Koopman learning can be used broadly for robust controller synthesis.

The existing algorithms assume a full state measurement and construct observables from that. With partial state observables like in biological systems, we construct a state space model for an output nonlinear difference equation model and identify a Koopman operator for that model. In Section II we briefly introduce the Koopman operator and DMD and the existing literature. In Section III we describe the experimental setup for obtaining the growth curve data of Pseudomonas putida by adding different concentrations of casein and glucose substrates to the media. In Section IV, we justify nonlinear autoregressive difference equation model is an appropriate choice for this system. In Section V, we formulate the Hankel DMD as a solution of the NAR model and bring out its issues and in Section VI, we formulate a state space model for the NAR model and use extended dynamic mode decomposition to identify a Koopman operator for the NAR model. In Section VII, we show that the algorithm is able to train a predictive Koopman operator, that predicts with 3.4% on the training data and 9% on the test data on extended forecasting tasks approximately 500 time steps ahead.

Consider a discrete-time autonomous nonlinear dynamical system

with $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is analytic. There exists a Koopman operator [34] of (1), which acts on a function space $\mathcal{F}$ as $\mathcal{K}$ : $\mathcal{F}$ $\rightarrow$ $\mathcal{F}$. This action can be given by

where the function $\psi:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is called an observable of the system and the set of all observables $\psi\triangleq\{\psi_{i}\}_{i=1}^{p},p\leq\infty$ on the system. Here $\mathcal{F}$ is invariant under the action of $\mathcal{K}$.

The most important property of the Koopman operator that we utilize is the linearity of the operator, in other words,

which follows from (2) since the composition operator is linear. Thus, we have that the Koopman operator of (1) is a linear operator that acts on observable functions $\psi(x_{k})$ and propagates them forward in time.

We describe the procedure adopted to obtain P. putida’s growth curve for varying concentrations of glucose and casein substrates in the media.

Data Collection: The microplate reader was set to $30^{o}$C and the shaker to 807 cycles per minute, with continuous double orbital mixing. The absorbance at 600 nanometers(nm), which is termed as the Optical Density at 600 nm ($OD_{600}$), was measured as a function of time for 27 hours. We assume in this work, as is widely accepted, that $OD_{600}$ measurements were collected in a linear regime, where cell population is proportional OD600 measurements. The obtained data along with the varying substrate concentrations are shown in Figure. 1.

The dynamics of the bacterial cell growth can be represented by

where the bacterial cell count ($N^{(b)}$), casein substrate concentration ($C$) and glucose substrate concentration ($G$) are the states of the system, $f$ is the nonlinear dynamics and $k$ is the discrete time index. We measure the $OD_{600}$ data as mentioned in section III and the output equation is given by

as $OD_{600}$ is a function of only the number of cells. Some of the existing empirical nonlinear models for growth curve dynamics include the Monod’s model [35] which uses a single substrate to form the foundation of the growth curve dynamics and in [36] and [37] multiple substrates are incorporated. Monod’s model is a two-state nonlinear dynamical system comprising of the substrate (S) and the number of bacteria ($N^{(b)}$):

where $r_{max}$ is the maximum growth rate and $K_{s}$ is the half velocity constant. As $N^{(b)}$ is the only variable of measurement in (5), we convert the model to a single differential equation containing only $N^{(b)}$

The existing models though heuristic, suggest that $N^{(b)}$ at any point in time is a function of the past

$N^{(b)}_{k}$ to be a function of its finite past. This is the general structure of the discrete nonlinear autoregressive (NAR) model given by

where the current output is a function of the past $\tau$ outputs.

Given the nonlinear system (4) with the state measurement given by (5) and modeled by the discrete time difference equation (8), Hankel DMD [38] is a suitable algorithm to solve the model identification problem with the NAR structure. The promising feature of using a DMD algorithm is that it identifies a linear state space representation which has a theoretical foundation in Koopman operator theory.

Given the autonomous state space system

where $x_{k}\in\mathbb{R}^{n}$ is the state, $\tilde{f}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is the dynamics, $y_{k}\in\mathbb{R}^{p}$ is the output and $h:\mathbb{R}^{n}\rightarrow\mathbb{R}^{p}$ is a nonlinear function that maps the state directly to itself, i.e. $x$ is identical to the output $y$, Hankel DMD constructs a Koopman model of the form

such that $\psi:\mathbb{R}^{p}\rightarrow\mathbb{R}^{N_{p}}$ is the dictionary of state inclusive observables of the state $\tilde{x}_{k}$ constructed by a nonlinear transformation of the corresponding output $y_{k}$ and $K$ is the Koopman operator. Regardless of full-state measurements, we nonetheless cast Hankel DMD in this form to compare it with our subsequent causal jump DMD algorithm.

Given the output measurements $\{y_{1},y_{2},..,y_{N}\}$, to identify an approximate Koopman operator $K$ using Hankel DMD, the time shifted Hankel matrices are constructed as

and the optimization problem

is solved using the Moore-Penrose pseudoinverse method mentioned in section II. This yields a solution of the form

Other than the last $N_{p}$ equations, the others are trivial. To construct an output predictor, we take the component $y_{k+\tau}$ of $\psi(y_{k+\tau})$ to get

where $\tilde{k}_{i}$ are the components of $k_{i}$ that map $\psi(y_{k+\tau-1})$ to $y_{k+\tau}$. More generally, this yields a nonlinear equation of the form

where the functions $\tilde{f}_{1},\tilde{f}_{2},\cdots,\tilde{f}_{\tau}$ have the same basis functions with different coefficients. This identifies a constrained NAR model as it imposes an additive structure on the basis of nonlinear models across time.

To identify a Koopman operator for the unconstrained NAR model (8), we formulate a state space representation for the NAR model with full state observation and identify an approximate Koopman operator for that model using the general class of DMD algorithms like extended DMD and deep DMD.

In this methodology, the problem is broken into two pieces —— the system identification aspect where we select the model structure and the dynamic mode decomposition aspect where we have to construct the dictionary of observables. We define a window parameter $\tau\in\mathbb{Z}_{>0}$ indicating how many past output snapshots are used to define a new extended dictionary of monomial observable functions, up to order $n_{o}\in\mathbb{Z}_{>0}$. The new $\tau$-dictionary defines a general extended dynamic mode decomposition problem, which we then solve using classical methods.

We proceed as follows: given the NAR model (8) with the system identification parameter $\tau$, we construct a state defined by

with $z_{k}\in\mathbb{R}^{p\tau}$.This yields the state space representation

where $\tilde{F}:\mathbb{R}^{p\tau}\rightarrow\mathbb{R}^{p\tau}$ represents the dynamics of the lifted ”state” model. The approximate EDMD model for the full output observable model is given by

where $\psi(z_{k})$ is the state inclusive dictionary of observables defined as

with $\varphi:\mathbb{R}^{p\tau}\rightarrow\mathbb{R}^{N_{p}}$ being a nonlinear transformation that constructs the nonlinear observables. Since the only additional information in the state $z_{k+1}$ when compared to the state $z_{k}$ is $y_{k+\tau+1}$, the output predictor form for the Koopman model can be identified considering the complete Koopman model and extracting the equation that corresponds to $y_{k+\tau+1}$ given by

The output predictor form keeps the general structure of the NAR model intact as opposed to the predictor identified by Hankel DMD which identified a constrained model. But, the issue with this model is the causality. It can be seen from (VI) that the Koopman model is non causal due to the overlap of outputs $y_{k}$ between the states $z_{k+1}$ and $z_{k}$. This identifies models that use future outputs to predict past outputs which are inadmissible as our system is causal. To identify a causal model, the property of (VI) proved in proposition 1 is very important.

To identify a causal Koopman model for the NAR system (8), we propagate the model (VI) by $\tau$ time steps to ensure no intersection of outputs between the states $z_{k+1}$ and $z_{k}$. We define a new state $x_{k}=z_{k\tau}$ which yields

Using proposition 1, we can say that the nonlinear state space model contains functions that are 1-step, 2-step, …, $\tau$-step predictors in the following form

where $f^{(i)}$ is the $i$-step predictor of the NAR model. We prove the existence of a Koopman operator for this model in Proposition 2.

Since the existence of a Koopman operator has been proved for the model (VI) in Proposition 2, we construct a state inclusive dictionary of observables

with $\varphi:\mathbb{R}^{p\tau}\rightarrow\mathbb{R}^{N_{p}}$ to define a Koopman model

This Koopman model is causal since there is no intersection of outputs between $x_{k+1}$ and $x_{k}$. The added feature of this model is that the DMD algorithm while identifying a Koopman operator, also simultaneously minimizes the 1-step, 2-step, …, $\tau$-step prediction error of the NAR model.

Now that we have a theoretical state space representation of a NAR model and established the conditions under which a Koopman operator exists, we turn our attention to the algorithm for identification of the Koopman operator. Given the data with M data sets and N data points in each data set $\{y^{(i)}_{1},y^{(i)}_{2},...,y^{(i)}_{N}\}$ where $i\in\{1,2,...M\}$ is the index of the data set, we construct the Hankel states $z_{k}$ and the dictionary of observables allowing the intermixing of states. We compile the observables into snapshot matrices $\tilde{\Psi}_{f}(z)$ and $\tilde{\Psi}_{p}(z)$ with a $\tau$ time step jump and solve the Koopman learning problem

using the methodology in Algorithm 1.

From the data-sets obtained in the plate reader experiments shown in Fig. 1, we used Algorithm 1 to implement extended DMD using monomials as the dictionary of observables

to identify an approximate Koopman operator for the state space model (VI) as a solution to the identification of the NAR model (8).

We use all the datasets in Fig. 1 to find a Koopman operator invariant to the substrate concentrations. They are broken equally into training, validation and test set. Given the two parameters $\tau$ (NAR model parameter) and $n_{o}$ (extended DMD parameter), we can find the optimal approximate Koopman operator by cumulatively iterating through the principal components and evaluating the summation of the mean squared error(MSE) of training and validation data. The number of principal components corresponding to the minimum MSE yields the optimal Koopman operator for a given $\tau$ and $n_{o}$. We then iterate through the two parameters to find the optimal model that minimizes the

By choosing $\tau=9$ and keeping the maximum order of monomials to 3, the Koopman operator has been identified and the prediction on the training data is shown in Figure 2 and it has an MSE of 3.4%. The identified Koopman operator has an MSE of 9% and the fit is shown in Figure 3.

The results on the experimental data suggest that Causal Jump DMD is a suitable candidate algorithm for identifying the Koopman operator of the population growth dynamics of bacteria and can also be extended in general to identify Koopman operators for NAR models.

In this paper, we introduced the microbial growth curve dynamics to motivate the usage of DMD algorithms to identify Koopman operators for NAR models. We formulated Hankel DMD as a state space representations of the NAR model and showed that it is restrictive in its structure. We construct a causal state space model for the NAR model and identify a Koopman operator for it using extended dynamic mode decomposition with a monomial dictionary of observables. We showed that it does a good job in predicting the population growth dynamics of Pseudomonas putida invariant to substrate concentrations. The future goals of this work is to use this model to identify the optimal media conditions for maximal and minimal growth of the microbe thereby enabling us to develop a general methodology to develop an external growth harness for microbes for dynamic growth control. To achieve this, we need to extend the mathematical models to allow for inputs and extend the identification to NARX and NARMAX models. Further, if we integrate this framework with deepDMD which aids in finding the observable functions in a parsimonious fashion, it renders a useful tool for identifying high dimensional linear models for nonlinear systems.

The authors gratefully acknowledge the funding of DARPA grants FA8750-17-C-0229, HR001117C0092, HR001117C0094, DEAC0576RL01830. The authors would also like to thank Professors Arun K. Tangirala, Igor Mezic, Milan Korda, Alexandre Mauroy, Nathan Kutz, Steve Haase, John Harer, Devin Strickland, and Eric Klavins for insightful discussions. Any opinions, findings, conclusions, or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the Defense Advanced Research Project Agency, the Department of Defense, or the United States government. This material is based on work supported by DARPA and AFRL under contract numbers FA8750-17-C-0229, HR001117C0092, HR001117C0094, DEAC0576RL01830.