Biomarkers enable objective monitoring of a given cell or state in a biological system and are widely used in research, biomanufacturing, and clinical practice.

However, identifying appropriate biomarkers that are both robustly measurable and capture a state accurately remains challenging.

We present a framework for biomarker identification based upon observability guided sensor selection.

Our methods, Dynamic Sensor Selection (DSS) and Structure-Guided Sensor Selection (SGSS), utilize temporal models and experimental data, offering a template for applying observability theory to unconventional data obtained from biological systems.

Unlike conventional methods that assume well-known, fixed dynamics, DSS adaptively select biomarkers or sensors that maximize observability while accounting for the time-varying nature of biological systems.

Additionally, SGSS incorporates structural information and diverse data to identify sensors which are resilient against inaccuracies in our model of the underlying system.

We validate our approaches by performing estimation on high dimensional systems derived from temporal gene expression data from partial observations.

§ INTRODUCTION
The selection of biomarkers is, at least in principle, a classic problem of systems theory.

Like many engineered, industrial, and socioeconomic processes, a central objective of experimental science lies in minimizing the requisite measurements and data collected, while preserving our capacity to accurately estimate, detect, and forecast the state of a complex system.

Historically, biomarker identification for cancer and disease has relied upon domain knowledge of the biological system <cit.>; however, such an approach is limited to explaining known or characterized phenomena <cit.>.

As recent advances in real time sequencing usher in a new era in genomics <cit.>, model-based biomarker selection has the potential to uncover uncommon sensors and biomarkers directly from [unbiased] data.

Model-based biomarker identification selects sensors to maximize the observability of dynamic models of the a biological system.

A system is called observable when the measurements or data collected from sensors provide sufficient information to determine unmeasured states <cit.>, and yet, while observability is a classic problem of systems theory <cit.>, many challenges remain in applying input/output models or state space models typical of controls engineering to the study of biological systems <cit.>.

In contrast to many physical systems – such as the pendulum, where the position and velocity capture the state and the equations of motion are known – dynamics of biological systems, and often the correct state representation, remain unknown.

The high dimensional and low temporal resolution of data gathered in many biological experiments present a further challenge, as these data are not compatible with standard methods for identification and learning dynamics of complex systems <cit.>.

In spite of these challenges, many models to predict cell trajectories during differentiation, perturbation, and reprogramming have been proposed <cit.>.

Remarkably, the landmark cell reprogramming (controller) experiments of Weintraub <cit.> and Yamanaka <cit.> were based on characterizing key biomarker genes (observer) of a target cell type.  This approach exemplifies a classic principle of control theory in biological systems: the dual concepts of controllability and observability.

Nevertheless, our observability analysis of gene regulation supports the notion that steering and monitoring biological systems are in fact not equivalent problems.

To address these challenges, we introduce a framework for biomarker selection founded on dynamic models of gene regulation.

We present two templates for sensor selection: Dynamic Sensor Selection (DSS) and Structure-Guided Sensor Selection (SGSS).

We demonstrate their efficacy in identifying biomarkers that optimize the observability of dynamics on gene regulatory networks derived from time-series gene expression datasets.

Our focus lies on the linear time variant (LTV) state-space model

\begin{equation}\label{eq: dynamics}
\begin{split}
%    \begin{cases}
\xvt{t+1}&=\Avt{t}\xvt{t}\\
\yvt{t}  &=\Cvt{t}\xvt{t}.
%    \end{cases}
\end{split}
\end{equation}

Here, $\xvt{t}\in\R^n$ is the system state, representing, for instance, the expression of each gene as a vector; $\Avt{t}\in\R^{n\times n}$ signifies a state transition matrix, akin to a gene regulatory network; $\Cvt{t}\in\R^{p_t\times n}$ stands as the sensor or measurement matrix, dictating our data collection process, so that $\yvt{t}\in\R^{p(t)}$ denotes our measurements or data (where $p(t)\ll n$).

If $\Av(t)$ is fixed for all $t$, we call the system linear time invariant (LTI).

When it is cost-prohibitive to measure the full state at each time, the sensor selection problem involves crafting a measurement matrix ($\Cvt{t}$) to ensure that the low-dimensional data ($\yvt{t}$) gathered throughout time or during an experiment offers the greatest insight into the complete state of the system ($\xvt{t}$).

at (0,0)

Mathematical Motivation.

The trajectory of three coupled Van der Pol oscillators is shown in 3D with each of the possible pairwise projections, obtained by observing only two oscillators, shown in a different plane. As the network shifts from random initial conditions, to synchronization and then phase locking, the amount of information gathered in each plane and utility of the different observations changes. The network model of oscillator coupling, synchronization dynamics, and governing equations of the system are shown at the top.

§ RESULTS

§.§ Dynamic Sensor Selection

Traditional methods for sensor selections first evaluate each variable as a sensor and then suggest monitoring as many top-ranked sensors as feasible, considering experimental constraints.

However, by alternating our sensors at each time step to measure different variables, the estimation of nonlinear and time-varying systems can be improved.

Motivating Example: Biological Oscillators.

As a first example, consider how to best observe a network of coupled oscillators

\begin{equation}\label{eq: oscillator network with observations}
\begin{split}
% \dot{\xvt{t}}
\dfrac{d\xvt{t}}{dt}& = \textbf{F}(\xvt{t}, \mu(t)) - \textbf{L}(\xvt{t})\\
\yvt{t} & = \Cvt{t}\xvt{t}.
\end{split}
\end{equation}

Here, $\xv$ is a vector representing the locations or values of each oscillator, $\Fv$ is the dynamics of individual oscillators with internal parameters $\mu,$ and $\Lv$ is the diffusion operator specific to the network structure.

Turing's theory of morphogenesis, Smale's two cell system, the repressilator, and other higher order motifs exemplify the dynamics of many biological systems described by <ref>, highlighting the importance of its observation <cit.>.

In <ref>, the trajectory of three coupled Van der Pol oscillators is shown.

When measuring the state of any two oscillators $\xv_1,\ \xv_2,$ or $\xv_3,$ the observations are 2D projections of the 3D trajectory.

With fixed sensors, the question “which oscillators are the best sensors?" is akin to asking “data from with 2D plane enables the best reconstruction of the 3D shape?"

This supposes the observed data measures the same two variables at all times.

However, the information content of each projection changes as the oscillators synchronize and phase lock (<ref>).

As a result, alternating the plane of observed data throughout time provides a clearer picture of the 3D shape and enables better estimation and prediction of the network trajectory.

As a root node (<ref>), $\xv_1$ is a good sensor to monitor the long term behavior of the system <cit.>.

Nevertheless, modification of oscillator connectivity or parametrization before reaching the limiting behavior may necessitate sensor reallocation.

As dynamics evolve, the number and distribution of sensors should change as well.

For example, synchronized networks require fewer sensors than unsynchronized ones (<ref>).

Similarly, changes in parameterization and connectivity of oscillators necessitates reallocation of sensors (SI.<ref>).

The cell cycle and differentiation stages exemplify temporal interactions where sensors are dynamically allocated.

For instance, interactions between key regulatory genes, such as P27, P21, CYCLIN D1, CDK4, and MYOD, change between proliferation, differentiation, and quiessence <cit.>.

The PIP-FUCCI biomarker, developed first as FUCCI, employs fluorescent biomarkers to distinguish cell cycle stages.

Initially, CDT1 and GEM gene expression distinguished the G1 stage from S, G2, and M <cit.> (SI.<ref>).

Adding PIP to monitor the PCNA gene enabled accurate detection of G2 phase transitions <cit.>.

Monitoring of CDT1, GEM, and PCNA between different cell cycle stages exemplifies DSS.

Recently, the introduction of adaptive sequencing, which allows for a sequencer to update in real time which genes, cells, or other markers are measured, provides a flexible framework for DSS on high dimensional genomics experiments <cit.>.

Maximizing Observability.

We propose two formulations of DSS.

Output energy measures the magnitude of the observation $\yvt{t}$ over time.

At time $T$, sensor selection to maximize energy $\mathcal{E}$ is formulated

\begin{equation}\label{eq: energy problem}
\max_{\Cvt{t}}\mathcal{E}\text{ for all $t,$ where }\mathcal{E}=\sum_{t=0}^T\yvt{t}^\top\yvt{t}.
\end{equation}

Adapting the approach of <cit.>, <ref> is solved through its Lagrange dual form (SI.<ref>).

<Ref> is predicated on the prediction of $\yvt{t},$ and while this assumption is reasonable in many scenarios

the observability Gramian offers a generalized measure of output energy.

To form the discrete-time observability Grammian, let $\Phi(t_2,t_1)= \Av_{t_2}\Av_{t_2-1}\cdots \Av_{t_1}$ denote the transition matrix from time $t_1$ to $t_2,$ so that the observability Gramian is

\begin{equation}\label{eq:gram}
\G_o=\sum_{t=0}^T\Phi(t,0)^\top\Cvt{t}^\top\Cvt{t}\Phi(t,0).
\end{equation}

Using the relation $\yvt{t}=\Cvt{t}\Phi(0,t)\xv_0,$ utilizing $\Cvt{t}\Phi(0,t)$ relaxes the need for a prediction of $\yvt{t}.$

By summing over the inner product of $\Cvt{t}\Phi(0,t)$ times its transpose, <ref> is a direct generalization of the energy $\mathcal{E}$ in <ref>.

In contrast to $\mathcal{E},$ $\G_o$ is a matrix rather than a scalar, and several measures of observability derived from $\G_o$ have been proposed. We consider the problem

\begin{equation}\label{eq: gramian optimization}
\min_{\Cvt{t}} J(\G_o)
\end{equation}

where $J(\cdot)$ denotes the trace, logarithm of the determinant, smallest eigenvalue, or rank, each of which provide a different observability measure (SI.<ref>).

For the trace, <ref> is solved with a linear program and can applied to high dimensional systems (SI.<ref>).

The methodologies of <ref> and <ref> can handle time-varying sensors, incorporate additional constraints such as SGSS, and support the implementation of scalable algorithms.

When compared with alternative sensor selection techniques in <ref>, these approaches are versatile.

Comparison of Methods.

Method     Criteria     DSS     SGSS     Targeted     Cost

Gramian     $\R$                 $\mathcal{O}(n^{2.5})$

Energy     $\R$                 $\mathcal{O}(kn^2)$

Structure     0/1                 $\mathcal{O}(e\sqrt{n})$

*$^*$ $n$ is the number of variables, $e$ is the number of interactions, and $k$ is the number of iterations in an eigenvalue calculation.

§.§ Structure Guided Sensor Selection

Sullivan's maxim “form ever proceeds function”
has long established the essence of the structure-function (S-F) causality dilemma.

When the system identification problem remains unresolved and the model of function contains errors, SGSS can exploit information in both the structure and function domains to constrain the DSS optimization problems.

SGSS considers system geometry and spatial arrangement, leveraging orthogonal experimental methods, to mitigate modeling errors and identify robust sensors.

Knowledge of the structure can aid our estimation and understanding of the dynamics based upon its function.

This perspective resonates with approaches in other domains both (1) algorithmically, where methods such as PageRank <cit.> and the the fast multipole method <cit.> leverage additional structures to compute on complex systems, and (2) from data, where the S-F relationship has been recognized in the brain <cit.>, gene regulation <cit.>, and community structures <cit.>.

Paired data of the position (S) and effect (F) together, such as genome structure (Hi-C, S) and gene expression (RNAseq, F), is more powerful than either information alone.

The observability can be viewed as either a binary or scalar feature.

A well known limitation of the binary Kalman observability test is that all systems in the form of <ref> are nearly observable.

Apropos of this constraint, in 1974, Lin proposed structural observability, where the sparsity structure of the operators $\Av$ and $\Cv$ determine a binary observability condition <cit.>.

In contrast, scalar measures of observability, derived directly from $\Avt{t}$ and $\Cvt{t},$ provided graded measures of observability (SI.<ref>).

While DSS adopts the scalar metric perspective, SGSS departs from the binary view of observability.

The structure considered by SGSS is independent of $\Avt{t}$'s sparsity but rather based upon external attributes or structures of our system that may not appear in the dynamics.

While this notion of structure in SGSS varies from Lin's usage of the word, the challenge remains the same: despite great experimental advancements over the past half century, system identification and learning the dynamics is not a solved problem for biological systems.

Obtaining the data for traditional system identification techniques to be successful is both experimentally challenging and cost-prohibitive.

Present methodologies have utilized LTI methods on time-series gene expression signals <cit.>, and SGSS seeks to complement these methods by incorporating readily accessible data pertaining to genome structure.

Experimental Framework. (A) Determining mathematical observability and performing sensor selection is typically considered within the gray box, where various optimization problems can be formulated based on the dynamics. In the broader experimental framework, we can observe two types of data: structural or functional. From the functional, time-series data, we learn dynamics. From the structural data, we can constrain our optimization. Together, DSS and SGSS should select sensors that are useful to estimate the internal system state. (B) Two systems with structural and functional relationships. B1: a small world network is shown where sensors (red) are placed evenly in space; B2: a simple network model (lower left triangle) is fit to Hi-C (upper left triangle) for chromosome 1. (C) Estimation results from the Proliferation and SWB25 datasets.

Observability in a Small World.

The tendency to meet strangers with mutual acquaintances is a byproduct of the spatial structures that shape small world networks.

For instance, Milgram's infamous experiment was guided by the geography of individuals from Nebraska to Boston <cit.>; the Watts-Strogatz (WS) model positions each vertex in a lattice before forming the network <cit.>; and the small world structure of gene regulatory networks is guided by the 3D organization of chromatin <cit.>.

In each case, the structure guides the formation and dissolution of interactions in the system.

The positioning of nodes on the lattice determines the expected value of each node as a sensor in small world generated with the WS model.

We constructed an ensemble of small world networks and evaluated the contribution of each node to the network observability based upon the Gramian.

The node contributions to observability on the lattices resembled their average contribution as sensors over all small world networks generated from each lattice (SI.<ref>).

Moreover, evenly spacing sensor nodes across the lattice proves an effective strategy for placing sensor nodes on small world networks when the precise small world adjacency structure $\Avt{t}$ is unknown (SI.<ref>).

This suggest that when the precise set or regulatory interactions or network edges of $\Avt{t}$ is only partially known but the underlying structure is well characterized sensor selection can be guided by the structure.

The Nucleus is a Small World.

While network models of gene regulation and chromatin architecture have been developed from self organization principles <cit.> and molecular dynamics simulations <cit.>, quantification of properties of the genome from structural data remains unexplored.

We developed a four parameter network model whose adjacency structure qualitatively mirrors Hi-C (<ref>, SI.<ref>).

and caveman properties capture the diagonal dominance and block structure characteristic of the fractal globule chromatin architecture and Hi-C data <cit.>.

Based upon our ability to fit networks to Hi-C with relatively few parameters, we proceeded to quantify the Small World Quotient (SWQ) for several Hi-C datasets.

Varying from individual chromosomes to the full genome, the SWQ of Hi-C networks was estimated at several different resolutions, and we observed small world properties in all cases.

The SWQ increased with the resolution and size of the Hi-C network and matrix, and the small world properties at multiple resolutions are consistent with the self similar or fractal structures are consistent with classic, multi-scale perspectives of Hi-C <cit.>.

Utilizing Hi-C data collected in parallel with the proliferation and reprogramming datasets, we evaluated the SWQ throughout time; however, in neither dataset did we observe a significant change in the SWQ throughout time.

The consistent small world propensity of Hi-C motivates augmenting DSS of gene regulatory networks based upon chromatin structure.

§.§ Application to Data

We applied DSS and SGSS on a range of data including both genomic and EEG signals (<ref>, SI.<ref>).

High dimensional, low frequency gene expression data are at the frontier of observability theory whereas the low dimensional, high frequency EEG signals are a classic problem to study.

We used standard approaches to learn LTI and LTV models (SI.<ref>), and the sensors of each model are assessed based on their ability to estimate the full system state from the sensor measurements (SI.<ref>).

Time series datasets.

Dataset   Dimension   Time Points   Replicates    Ref.

Proliferation     19235      8     2   <cit.>

Reprogramming   19235    15     3    <cit.>

myogenicSignal       404     15     3

SBW25       624      9     2   <cit.>

EEG   64   160   109   <cit.>

To validate our models of gene expression dynamics, we employed established biomarkers from the literature to estimate gene expression during cell proliferation <cit.> (SI.<ref>).

Human fibroblasts were synchronized in terms of both the cell cycle stage and circadian rhythm, offering optimal conditions for learning LTI and LTV models.

For sensor selection, we employed the KEGG pathway database, which contains manually curated sets of genes <cit.> (SI.<ref>).

Initially, we investigate pathways associated with the cell cycle, such as the Basal Transcription Factors (hsa03022), Cell Cycle (hsa04110), Circadian Rhythm (hsa04710), Circadian Entrainment (hsa04713), and Cellular Senescence (hsa04218) pathways.

LTV models had median component wise errors bounded near 10%, which outperformed LTI models when using sensors from all pathways except hsa04713 (SI.fig.<ref>).

Although LTV dynamics generally offer superior estimation, we observed that they exhibit decreased robustness due issues such as over fitting and poor conditioning of the observability matrix (SI.<ref>).

Considering the role of transcription factors (TFs) in determining cell fate and the duality of controllability and observability, we hypothesized that including TFs is essential to forming effective sensor sets.

Consistent with this, while hsa04713 contained the third most genes of the sensor pathways considered thus far, it contained no TFs.

Repeating the estimation anew with all human KEGG pathways as sensor sets ($n=346$), we discovered neither the presence of a large number of TFs nor a large sensor set are necessary for good estimation, thereby challenging our hypothesis (<ref>, SI.<ref>).

Mathematically, TFs' effectiveness as controllers but not observers, which is contrary to linear systems theory, is ascribed to the nonlinearity of biological systems.

Biologically, TFs' relatively low expression levels result in low output energy and less variability in their concentrations, necessitating more sensitive observer and estimation approaches.

We observed a bifurcating behaviour in the estimation procedure.

Of the sensors that poorly estimate the initial state, the failed predictions deviate from biologically meaningful values by several orders of magnitude.

This improves the interpretability of our approach by offering a clear indicator of failure, even in cases where the true state of the system is unknown.

Pesticide Detection.

We build models of the gene regulatory network for Pseudomonas fluorescens SBW25 and selected biomarkers for malathion detection, a commonly used insecticide <cit.>.

In one model, we learned LTI dynamics ($\Av$) and time invariant sensors ($\Cv$), and in another model, we learned LTV dynamics ($\Avt{t}$) and used DSS ($\Cvt{t}$).

Varying the number of sensors, we assessed the estimation capabilities of each model, and found that LTV dynamics and DSS improved prediction accuracy for reconstructing the expression levels of individual genes.

Although <ref> and <ref> can always be further maximized by adding more sensors, in practice, increasing the number of sensors may not improve estimation, as illustrated in <ref>.C1.

Cellular Reprogramming.

The low efficiency of Weintraub's famous myogenic reprogramming experiment remains an active challenge in cell reprogramming <cit.> (SI.<ref>).

Monitoring cells throughout reprogramming may offer insight to this issue;

however, both formulations of DSS fail to perform well on this system, likely due to the asynchronized and noisy experimental conditions.

We applied SGSS to improve state estimation and increase observability by selecting spatially distributed genes.

Based on the hypothesis that colocalized genes are coregulated, we clustered genes according to Hi-C data and constrained DSS to select at most one sensor from each cluster (SI.<ref>).

By including constrained selection from Hi-C, the distribution of sensors across chromosomes shifted to mirror the distribution of genes (SI.<ref>-<ref>).

While we cannot measure the spatial proximity of clustered genes, we observed correlation in the expression values of several gene clusters, consistent with the concept of transcription factories.

Regardless, the estimation was improved by the Hi-C constrained SGSS.

When using few sensors, SGSS reduced the variance and improved estimation accuracy by approximately 25%.

To improve estimation further, we amplified the weak reprogramming signal by sampling genes involved in myogenesis and proliferation (SI.<ref>).

This targeted dataset provides improved conditions for biomarker identification, counteracting the experimental conditions of reprogramming.

Under these conditions, the estimation of the initial state for the reduced data shows median component-wise errors below 15% with all combinations of fixed or dynamic sensors from energy or Gramian based selection.

To close the design $\xrightarrow{}$ build $\xrightarrow{}$ test loop for myogenicSignal, we utilized these sensors to estimate the state of the complete reprogramming data.

Sensors selected from the reduced data, when optimized for energy, fail to estimate the full data well.

This occurs since the high energy genes in the targeted data have low energy in the complete reprogramming signals.

However, biomarkers identified via the Gramian on the targeted dataset continue to perform well at estimating the full data.

The median component-wise error is improved when applying Gramian selected biomarkers from the targeted data to the full data.

Converse of targeted observability <cit.>, where sensors are selected on the full reprogramming time series to observe only the myogenic signal, Gramian based sensor selection identified genes on the reduced data that estimate the full system well.

Beyond Genomics.
We employed DSS to rank different sensors observed in EEG signals.

The brain's properties are well-documented, and current research suggest EEGs are observabile with few sensors <cit.>.

We ranked the sensors of 64-lead EEG signals based on their contributions to output energy and the observability Gramian.

Relative to the genomics data, where synchronized or controled experiments have low frequency, high dimensional measurements, EEG data are high frequency and low dimensional, and the EEG signals are unsynchronized.

Instead EEG signals were partitioned according to different tasks the participants performed, such as opening or closing their eyes, prior to performing sensor selection.

The sensor ranking exhibit great variability across different activities, which underscores the utility of DSS, when participants change between tasks, a common occurrence in clinical settings.

In this context, the significance of sensors is determined by the participants' activities or states rather than specific time points from the start of the EEG signals.

Consistent with the principles of DSS, transitions between states coincide with variations in the most relevant sensors.

§ DISCUSSION

Many biological systems exhibit high dimensional, unknown dynamics that evolve overtime, often in an unpredictable manner.

Here, we have extended state space and network observability methods to develop a template for the observability of systems that are constrained to omit high dimensional and low temporally resolved data.

Beyond the initial step of over measuring the system prior to sensor selection, we stress the assumptions and limitations of our study.

In particular, time dependent observability is sensible when monitoring  synchronized or perturbed systems, where control signals act as reference points in time for sensor selection.

When dealing with systems where the state evolves but cannot be determined a priori based upon the time, it is more appropriate to consider state dependent observability.

Furthermore, while we apply SGSS based upon gene clustering from Hi-C data to identify transcription factors, alternative procedures based upon gene regulatory networks, chromatin accessibility, or alternative data and clustering techniques may be utilized.

Although our application of SGSS based on Hi-C improves estimation, several user defined choices are made in this process which can be further refined.

The inherent flexibility and freedom of these procedures to be adapted for different systems and data make these templates versatile for sensor selection both in and beyond the genome.

Our work also raises several questions worthy of future pursuit.

Implicit in our state space model is the representation of genes as model states.

Expanding the state space representation to incorporate isoform, chromatin accessibility, or other exponentially available omics data could enhance these models.

Moreover, the time series experimental datasets utilized in our study are divorced from RNA veolocity or pseudotime approaches that are also utilized to study genome dynamics.

Such methodologies may be married with the framework of our study to facilitate the analysis of single cell resolved dynamics.

As contemporary trends in both science and industry emphasize harnessing computing power for modeling from larger data, it's crucial to highlight that data quantity must not compromise focused experimentation.

More data, Big Data, and recent excitement around AI models are not a panacea for science.

Rather the collect of data to maximize observability must work in parsimony with modeling approaches to gain new insights to complex systems.

§ ACKNOWLEDGMENTS
We would like to thank the members of the Rajapakse Lab for helpful and inspiring discussion. This work is supported by Air Force Office of Scientific Research (AFOSR) under award number FA9550-22-1-0215 (IR) and NIGMS GM150581 (JP).

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleEmergence of dynamic properties in network hypermotifs

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleEndocytosis as a stabilizing mechanism for tissue homeostasis

journaltitleSIAM Journal on Applied Mathematics
titleBifurcation analysis of an SIRS epidemic model with generalized incidence

titleAn introduction to systems biology: design principles of biological circuits

Nature Publishing Group UK London

journaltitleNature Reviews Genetics
titleNetwork motifs: theory and experimental approaches

Chalmers Tekniska Hogskola (Sweden)

titleNonlinear Observability and Identi ability: General Theory and a Case Study of a Kinetic Model for S. cerevisiae

journaltitleIEEE Transactions on Network Science and Engineering
titleStructural accessibility and structural observability of nonlinear networked systems

titleSystem identification—a survey

journaltitleNonlinear Analysis: Theory, Methods & Applications
titleControllability and observability of polynomial dynamical systems

Sage Publications Sage CA: Los Angeles, CA

journaltitleThe Neuroscientist
titleSmall-world brain networks revisited

Sage Publications Sage CA: Thousand Oaks, CA

journaltitleThe neuroscientist
titleSmall-world brain networks

Nature Publishing Group UK London

journaltitleNature communications
titleComplex small-world regulatory networks emerge from the 3D organisation of the human genome

National Acad Sciences

journaltitleProceedings of the national academy of sciences
titleDiscovering governing equations from data by sparse identification of nonlinear dynamical systems

booktitle2013 American control conference
titleOn strong structural controllability of networked systems: A constrained matching approach

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleFunctional organization of the human 4D Nucleome

Nature Publishing Group UK London

titleLive-seq enables temporal transcriptomic recording of single cells

Nature Publishing Group UK London

journaltitleScientific reports
titlePolymer physics of chromosome large-scale 3D organisation

Public Library of Science San Francisco, USA

journaltitlePloS one
titleNodal dynamics, not degree distributions, determine the structural controllability of complex networks

titleOptimal estimation of dynamic systems

Cold Spring Harbor Laboratory

titleAutomated In Vitro Wound Healing Assay

Princeton University Press Princeton, NJ

titleBiomolecular feedback systems

Nature Publishing Group UK London

journaltitleNature Communications
titleDeciphering multi-way interactions in the human genome

Nature Publishing Group UK London

titleA synthetic oscillatory network of transcriptional regulators

journaltitlePubl. math. debrecen
titleOn random graphs I

journaltitlePubl. math. inst. hung. acad. sci
titleOn the evolution of random graphs

Nature Publishing Group UK London

journaltitleNature reviews cancer
titleThe case for early detection

SAGE Publications Sage UK: London, England

titleSmall world sampling: Qualitative sample reliability and validity for efficient and effective recruitment of journalists as research participants

journaltitleBiophysical journal
titleImpulses and physiological states in theoretical models of nerve membrane

Oxford University Press

titlePanglaoDB: a web server for exploration of mouse and human single-cell RNA sequencing data

Nature Publishing Group UK London

titleConstruction of a genetic toggle switch in Escherichia coli

Oldenbourg Wissenschaftsverlag

titleOn global and local observability of nonlinear polynomial systems: A decidable criterion

journaltitleThe Annals of Mathematical Statistics
titleRandom graphs

journaltitleJournal of the Society for Industrial and Applied Mathematics, Series A: Control
titleControllability and observability in multivariable control systems

National Acad Sciences

journaltitleProceedings of the national academy of sciences
titleCommunity structure in social and biological networks

journaltitlesiam REVIEW
titlePageRank beyond the web

journaltitleActa Biotheoretica
titleThe Goodwin oscillator and its legacy

journaltitleAdvances in enzyme regulation
titleOscillatory behavior in enzymatic control processes

titleTemporal organization in cells. A dynamic theory of cellular control processes.

journaltitleSociological theory
titleThe strength of weak ties: A network theory revisited

University of Chicago Press

journaltitleAmerican journal of sociology
titleThe strength of weak ties

Taylor & Francis

journaltitleCell Cycle
titleAccurate delineation of cell cycle phase transitions in living cells with PIP-FUCCI

journaltitleIEEE Transactions on Circuits and Systems
titleDynamics of the van der Pol equation

journaltitleSIAM Journal on Applied Dynamical Systems
titleChaos in the Hodgkin–Huxley model

American Association for the Advancement of Science

titleCorrelation of terminal cell cycle arrest of skeletal muscle with induction of p21 by MyoD

Nature Publishing Group UK London

journaltitleNature biotechnology
titleCancer biomarkers: a systems approach

Nature Publishing Group UK London

numberSuppl 6761
titleFrom molecular to modular cell biology

Nature Publishing Group UK London

journaltitleNature Communications
titleLearning perturbation-inducible cell states from observability analysis of transcriptome dynamics

CUP Archive

titleTheory and applications of Hopf bifurcation

journaltitleIEEE Transactions on automatic control
titleNonlinear controllability and observability

Elsevier New York, NY, USA

journaltitleEncyclopedia of mathematical physics
titleWeakly coupled oscillators

journaltitleJournal of molecular biology
titleGenetic regulatory mechanisms in the synthesis of proteins

Cold Spring Harbor Laboratory Press

booktitleCold Spring Harbor symposia on quantitative biology
titleOn the regulation of gene activity

journaltitleJournal of the Society for Industrial and Applied Mathematics, Series A: Control
titleMathematical description of linear dynamical systems

Oxford University Press

journaltitleNucleic acids research
titleKEGG: new perspectives on genomes, pathways, diseases and drugs

booktitleProceedings of the thirty-second annual ACM symposium on Theory of computing
titleThe small-world phenomenon: An algorithmic perspective

Sage Publications Sage CA: Los Angeles, CA

journaltitleAmerican sociological review
titleThe small world of Germany and the durability of national networks

American Institute of Physics

journaltitleThe Journal of Chemical Physics
titleDynamics of entangled linear polymer melts: A molecular-dynamics simulation

booktitleInternational Symposium on Mathematical Problems in Theoretical Physics: January 23–29, 1975, Kyoto University, Kyoto/Japan
titleSelf-entrainment of a population of coupled non-linear oscillators

titleDynamic mode decomposition: data-driven modeling of complex systems

titleThe human transcription factors

Wiley New York

titleFoundations of optimal control theory

Annual Reviews

journaltitleAnnual Review of Control, Robotics, and Autonomous Systems
titleFast and flexible multiagent decision-making

American Association for the Advancement of Science

titleThe fundamental advantages of temporal networks

Nature Publishing Group US New York

journaltitleNature Biotechnology
titleMulti-omic single-cell velocity models epigenome–transcriptome interactions and improves cell fate prediction

American Association for the Advancement of Science

titleComprehensive mapping of long-range interactions reveals folding principles of the human genome

journaltitleIEEE Transactions on Automatic Control
titleStructural controllability

titleGenome architecture mediates transcriptional control of human myogenic reprogramming

Nature Publishing Group UK London

titleControllability of complex networks

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleObservability of complex systems

Prentice Hall PTR

seriesPrentice Hall information and system sciences series
titleSystem Identification: Theory for the User

journaltitleAnnual Reviews in Control
titlePerspectives on system identification

booktitleSignal analysis and prediction
titleSystem identification

journaltitlearXiv preprint arXiv:2308.02861
titleA unified-field theory of genome organization and gene regulation

journaltitleThe Hopf bifurcation and its applications
titleA mathematical model of two cells via Turing’s equation

Springer Science & Business Media

titleThe Hopf bifurcation and its applications

journaltitleSIAM Journal on Control and Optimization
titleStrong structural controllability

New York

journaltitlePsychology today
titleThe small world problem

American Association for the Advancement of Science

titleNetwork motifs: simple building blocks of complex networks

journaltitlearXiv preprint arXiv:2401.16372
titleDuality between controllability and observability for target control and estimation in networks

journaltitleIEEE Control Systems Letters
titleTarget controllability and target observability of structured network systems

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleFunctional observability and target state estimation in large-scale networks

journaltitleIEEE Transactions on Automatic Control
titleOn the structural and strong structural controllability of undirected networks

booktitleHandbook of Materials Modeling: Methods
titleThe Bond Fluctuation Model and Other Lattice Models

journaltitleProceedings of the IRE
titleAn active pulse transmission line simulating nerve axon

Cambridge University Press

journaltitleNetwork Science
titleHow small is it? Comparing indices of small worldliness

journaltitleJournal of Statistical Physics
titleModels of the small world

National Acad Sciences

journaltitleProceedings of the national academy of sciences
titleModularity and community structure in networks

booktitleProceedings of the 2005, American Control Conference, 2005.
titleUltrafast consensus in small-world networks

journaltitleIEEE Transactions on Control of Network Systems
titleControllability metrics, limitations and algorithms for complex networks

booktitle2015 54th IEEE Conference on Decision and Control (CDC)
titleMinimum number of probes for brain dynamics observability

Gauthier-Villars et fils, imprimeurs-libraires

titleLes méthodes nouvelles de la mécanique céleste

titleHyperstability of control systems

titleMapping transcriptomic vector fields of single cells

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleDynamics and control of state-dependent networks for probing genomic organization

John Wiley & Sons, Ltd Chichester, UK

journaltitleMolecular systems biology
titleNetworking the nucleus

titleA 3D map of the human genome at kilobase resolution reveals principles of chromatin looping

journaltitleJournal of computational physics
titleRapid solution of integral equations of classical potential theory

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleAlgorithm for cellular reprogramming

titleVisualizing spatiotemporal dynamics of multicellular cell-cycle progression

journaltitleIEEE Transactions on biomedical engineering
titleBCI2000: a general-purpose brain-computer interface (BCI) system

Annual Reviews

journaltitleAnnual Review of Fluid Mechanics
titleDynamic mode decomposition and its variants

booktitleProceedings of the 2001 international symposium on Symbolic and algebraic computation
titleA probabilistic algorithm to test local algebraic observability in polynomial time

journaltitleIEEE Transactions on Automatic control
titleStructural controllability of multiinput linear systems

journaltitleSystems & Control Letters
titleA concept of local observability

journaltitleIEEE spectrum
titleLeast-squares estimation: from Gauss to Kalman

titleThe proliferation-quiescence decision is controlled by a bifurcation in CDK2 activity at mitotic exit

titleThe small world of the cerebral cortex

journaltitleIEEE Transactions on Control of Network Systems
titleOn submodularity and controllability in complex dynamical networks

titleInduction of pluripotent stem cells from adult human fibroblasts by defined factors

Mary Ann Liebert, Inc. 140 Huguenot Street, 3rd Floor New Rochelle, NY 10801 USA

journaltitleBRAIN CONNECTIVITY
titleThe ubiquity of small-world networks

journaltitleBulletin of mathematical biology
titleThe chemical basis of morphogenesis

Nature Publishing Group UK London

titleUltra-fast deep-learned CNS tumour classification during surgery

titleObservability and structural identifiability of nonlinear biological systems

WW Norton & Company

titleSix degrees: The science of a connected age

Nature Publishing Group

titleCollective dynamics of ‘small-world’networks

Nature Publishing Group US New York

journaltitleNature Biotechnology
titleDynamic, adaptive sampling during nanopore sequencing using Bayesian experimental design

titleThe MyoD family and myogenesis: redundancy, networks, and thresholds

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleActivation of muscle-specific genes in pigment, nerve, fat, liver, and fibroblast cell lines by forced expression of MyoD.

IOP Publishing

journaltitlePhysical biology
titleThe Goodwin model revisited: Hopf bifurcation, limit-cycle, and periodic entrainment

Nature Publishing Group UK London

titleElite and stochastic models for induced pluripotent stem cell generation

[University of New Mexico, University of Chicago Press]

abstractData from a voluntary association are used to construct a new formal model for a traditional anthropological problem, fission in small groups. The process leading to fission is viewed as an unequal flow of sentiments and information across the ties in a social network. This flow is unequal because it is uniquely constrained by the contextual range and sensitivity of each relationship in the network. The subsequent differential sharing of sentiments leads to the formation of subgroups with more internal stability than the group as a whole, and results in fission. The Ford-Fulkerson labeling algorithm allows an accurate prediction of membership in the subgroups and of the locus of the fission to be made from measurements of the potential for information flow across each edge in the network. Methods for measurement of potential information flow are discussed, and it is shown that all appropriate techniques will generate the same predictions.
journaltitleJournal of Anthropological Research
titleAn Information Flow Model for Conflict and Fission in Small Groups

Taylor & Francis

journaltitleDynamics and stability of systems
titleHopf bifurcations in competitive three-dimensional Lotka–Volterra systems

titleCircuit design features of a stable two-cell system

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleEmergence of dynamic properties in network hypermotifs

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleEndocytosis as a stabilizing mechanism for tissue homeostasis

journaltitleSIAM Journal on Applied Mathematics
titleBifurcation analysis of an SIRS epidemic model with generalized incidence

titleAn introduction to systems biology: design principles of biological circuits

Nature Publishing Group UK London

journaltitleNature Reviews Genetics
titleNetwork motifs: theory and experimental approaches

Chalmers Tekniska Hogskola (Sweden)

titleNonlinear Observability and Identi ability: General Theory and a Case Study of a Kinetic Model for S. cerevisiae

journaltitleIEEE Transactions on Network Science and Engineering
titleStructural accessibility and structural observability of nonlinear networked systems

titleSystem identification—a survey

journaltitleNonlinear Analysis: Theory, Methods & Applications
titleControllability and observability of polynomial dynamical systems

Sage Publications Sage CA: Los Angeles, CA

journaltitleThe Neuroscientist
titleSmall-world brain networks revisited

Sage Publications Sage CA: Thousand Oaks, CA

journaltitleThe neuroscientist
titleSmall-world brain networks

Nature Publishing Group UK London

journaltitleNature communications
titleComplex small-world regulatory networks emerge from the 3D organisation of the human genome

National Acad Sciences

journaltitleProceedings of the national academy of sciences
titleDiscovering governing equations from data by sparse identification of nonlinear dynamical systems

booktitle2013 American control conference
titleOn strong structural controllability of networked systems: A constrained matching approach

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleFunctional organization of the human 4D Nucleome

Nature Publishing Group UK London

titleLive-seq enables temporal transcriptomic recording of single cells

Nature Publishing Group UK London

journaltitleScientific reports
titlePolymer physics of chromosome large-scale 3D organisation

Public Library of Science San Francisco, USA

journaltitlePloS one
titleNodal dynamics, not degree distributions, determine the structural controllability of complex networks

titleOptimal estimation of dynamic systems

Cold Spring Harbor Laboratory

titleAutomated In Vitro Wound Healing Assay

Princeton University Press Princeton, NJ

titleBiomolecular feedback systems

Nature Publishing Group UK London

journaltitleNature Communications
titleDeciphering multi-way interactions in the human genome

Nature Publishing Group UK London

titleA synthetic oscillatory network of transcriptional regulators

journaltitlePubl. math. debrecen
titleOn random graphs I

journaltitlePubl. math. inst. hung. acad. sci
titleOn the evolution of random graphs

Nature Publishing Group UK London

journaltitleNature reviews cancer
titleThe case for early detection

SAGE Publications Sage UK: London, England

titleSmall world sampling: Qualitative sample reliability and validity for efficient and effective recruitment of journalists as research participants

journaltitleBiophysical journal
titleImpulses and physiological states in theoretical models of nerve membrane

Oxford University Press

titlePanglaoDB: a web server for exploration of mouse and human single-cell RNA sequencing data

Nature Publishing Group UK London

titleConstruction of a genetic toggle switch in Escherichia coli

Oldenbourg Wissenschaftsverlag

titleOn global and local observability of nonlinear polynomial systems: A decidable criterion

journaltitleThe Annals of Mathematical Statistics
titleRandom graphs

journaltitleJournal of the Society for Industrial and Applied Mathematics, Series A: Control
titleControllability and observability in multivariable control systems

National Acad Sciences

journaltitleProceedings of the national academy of sciences
titleCommunity structure in social and biological networks

journaltitlesiam REVIEW
titlePageRank beyond the web

journaltitleActa Biotheoretica
titleThe Goodwin oscillator and its legacy

journaltitleAdvances in enzyme regulation
titleOscillatory behavior in enzymatic control processes

titleTemporal organization in cells. A dynamic theory of cellular control processes.

journaltitleSociological theory
titleThe strength of weak ties: A network theory revisited

University of Chicago Press

journaltitleAmerican journal of sociology
titleThe strength of weak ties

Taylor & Francis

journaltitleCell Cycle
titleAccurate delineation of cell cycle phase transitions in living cells with PIP-FUCCI

journaltitleIEEE Transactions on Circuits and Systems
titleDynamics of the van der Pol equation

journaltitleSIAM Journal on Applied Dynamical Systems
titleChaos in the Hodgkin–Huxley model

American Association for the Advancement of Science

titleCorrelation of terminal cell cycle arrest of skeletal muscle with induction of p21 by MyoD

Nature Publishing Group UK London

journaltitleNature biotechnology
titleCancer biomarkers: a systems approach

Nature Publishing Group UK London

numberSuppl 6761
titleFrom molecular to modular cell biology

Nature Publishing Group UK London

journaltitleNature Communications
titleLearning perturbation-inducible cell states from observability analysis of transcriptome dynamics

CUP Archive

titleTheory and applications of Hopf bifurcation

journaltitleIEEE Transactions on automatic control
titleNonlinear controllability and observability

Elsevier New York, NY, USA

journaltitleEncyclopedia of mathematical physics
titleWeakly coupled oscillators

journaltitleJournal of molecular biology
titleGenetic regulatory mechanisms in the synthesis of proteins

Cold Spring Harbor Laboratory Press

booktitleCold Spring Harbor symposia on quantitative biology
titleOn the regulation of gene activity

journaltitleJournal of the Society for Industrial and Applied Mathematics, Series A: Control
titleMathematical description of linear dynamical systems

Oxford University Press

journaltitleNucleic acids research
titleKEGG: new perspectives on genomes, pathways, diseases and drugs

booktitleProceedings of the thirty-second annual ACM symposium on Theory of computing
titleThe small-world phenomenon: An algorithmic perspective

Sage Publications Sage CA: Los Angeles, CA

journaltitleAmerican sociological review
titleThe small world of Germany and the durability of national networks

American Institute of Physics

journaltitleThe Journal of Chemical Physics
titleDynamics of entangled linear polymer melts: A molecular-dynamics simulation

booktitleInternational Symposium on Mathematical Problems in Theoretical Physics: January 23–29, 1975, Kyoto University, Kyoto/Japan
titleSelf-entrainment of a population of coupled non-linear oscillators

titleDynamic mode decomposition: data-driven modeling of complex systems

titleThe human transcription factors

Wiley New York

titleFoundations of optimal control theory

Annual Reviews

journaltitleAnnual Review of Control, Robotics, and Autonomous Systems
titleFast and flexible multiagent decision-making

American Association for the Advancement of Science

titleThe fundamental advantages of temporal networks

Nature Publishing Group US New York

journaltitleNature Biotechnology
titleMulti-omic single-cell velocity models epigenome–transcriptome interactions and improves cell fate prediction

American Association for the Advancement of Science

titleComprehensive mapping of long-range interactions reveals folding principles of the human genome

journaltitleIEEE Transactions on Automatic Control
titleStructural controllability

titleGenome architecture mediates transcriptional control of human myogenic reprogramming

Nature Publishing Group UK London

titleControllability of complex networks

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleObservability of complex systems

Prentice Hall PTR

seriesPrentice Hall information and system sciences series
titleSystem Identification: Theory for the User

journaltitleAnnual Reviews in Control
titlePerspectives on system identification

booktitleSignal analysis and prediction
titleSystem identification

journaltitlearXiv preprint arXiv:2308.02861
titleA unified-field theory of genome organization and gene regulation

journaltitleThe Hopf bifurcation and its applications
titleA mathematical model of two cells via Turing’s equation

Springer Science & Business Media

titleThe Hopf bifurcation and its applications

journaltitleSIAM Journal on Control and Optimization
titleStrong structural controllability

New York

journaltitlePsychology today
titleThe small world problem

American Association for the Advancement of Science

titleNetwork motifs: simple building blocks of complex networks

journaltitlearXiv preprint arXiv:2401.16372
titleDuality between controllability and observability for target control and estimation in networks

journaltitleIEEE Control Systems Letters
titleTarget controllability and target observability of structured network systems

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleFunctional observability and target state estimation in large-scale networks

journaltitleIEEE Transactions on Automatic Control
titleOn the structural and strong structural controllability of undirected networks

booktitleHandbook of Materials Modeling: Methods
titleThe Bond Fluctuation Model and Other Lattice Models

journaltitleProceedings of the IRE
titleAn active pulse transmission line simulating nerve axon

Cambridge University Press

journaltitleNetwork Science
titleHow small is it? Comparing indices of small worldliness

journaltitleJournal of Statistical Physics
titleModels of the small world

National Acad Sciences

journaltitleProceedings of the national academy of sciences
titleModularity and community structure in networks

booktitleProceedings of the 2005, American Control Conference, 2005.
titleUltrafast consensus in small-world networks

journaltitleIEEE Transactions on Control of Network Systems
titleControllability metrics, limitations and algorithms for complex networks

booktitle2015 54th IEEE Conference on Decision and Control (CDC)
titleMinimum number of probes for brain dynamics observability

Gauthier-Villars et fils, imprimeurs-libraires

titleLes méthodes nouvelles de la mécanique céleste

titleHyperstability of control systems

titleMapping transcriptomic vector fields of single cells

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleDynamics and control of state-dependent networks for probing genomic organization

John Wiley & Sons, Ltd Chichester, UK

journaltitleMolecular systems biology
titleNetworking the nucleus

titleA 3D map of the human genome at kilobase resolution reveals principles of chromatin looping

journaltitleJournal of computational physics
titleRapid solution of integral equations of classical potential theory

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleAlgorithm for cellular reprogramming

titleVisualizing spatiotemporal dynamics of multicellular cell-cycle progression

journaltitleIEEE Transactions on biomedical engineering
titleBCI2000: a general-purpose brain-computer interface (BCI) system

Annual Reviews

journaltitleAnnual Review of Fluid Mechanics
titleDynamic mode decomposition and its variants

booktitleProceedings of the 2001 international symposium on Symbolic and algebraic computation
titleA probabilistic algorithm to test local algebraic observability in polynomial time

journaltitleIEEE Transactions on Automatic control
titleStructural controllability of multiinput linear systems

journaltitleSystems & Control Letters
titleA concept of local observability

journaltitleIEEE spectrum
titleLeast-squares estimation: from Gauss to Kalman

titleThe proliferation-quiescence decision is controlled by a bifurcation in CDK2 activity at mitotic exit

titleThe small world of the cerebral cortex

journaltitleIEEE Transactions on Control of Network Systems
titleOn submodularity and controllability in complex dynamical networks

titleInduction of pluripotent stem cells from adult human fibroblasts by defined factors

Mary Ann Liebert, Inc. 140 Huguenot Street, 3rd Floor New Rochelle, NY 10801 USA

journaltitleBRAIN CONNECTIVITY
titleThe ubiquity of small-world networks

journaltitleBulletin of mathematical biology
titleThe chemical basis of morphogenesis

Nature Publishing Group UK London

titleUltra-fast deep-learned CNS tumour classification during surgery

titleObservability and structural identifiability of nonlinear biological systems

WW Norton & Company

titleSix degrees: The science of a connected age

Nature Publishing Group

titleCollective dynamics of ‘small-world’networks

Nature Publishing Group US New York

journaltitleNature Biotechnology
titleDynamic, adaptive sampling during nanopore sequencing using Bayesian experimental design

titleThe MyoD family and myogenesis: redundancy, networks, and thresholds

National Acad Sciences

journaltitleProceedings of the National Academy of Sciences
titleActivation of muscle-specific genes in pigment, nerve, fat, liver, and fibroblast cell lines by forced expression of MyoD.

IOP Publishing

journaltitlePhysical biology
titleThe Goodwin model revisited: Hopf bifurcation, limit-cycle, and periodic entrainment

Nature Publishing Group UK London

titleElite and stochastic models for induced pluripotent stem cell generation

[University of New Mexico, University of Chicago Press]

abstractData from a voluntary association are used to construct a new formal model for a traditional anthropological problem, fission in small groups. The process leading to fission is viewed as an unequal flow of sentiments and information across the ties in a social network. This flow is unequal because it is uniquely constrained by the contextual range and sensitivity of each relationship in the network. The subsequent differential sharing of sentiments leads to the formation of subgroups with more internal stability than the group as a whole, and results in fission. The Ford-Fulkerson labeling algorithm allows an accurate prediction of membership in the subgroups and of the locus of the fission to be made from measurements of the potential for information flow across each edge in the network. Methods for measurement of potential information flow are discussed, and it is shown that all appropriate techniques will generate the same predictions.
journaltitleJournal of Anthropological Research
titleAn Information Flow Model for Conflict and Fission in Small Groups

Taylor & Francis

journaltitleDynamics and stability of systems
titleHopf bifurcations in competitive three-dimensional Lotka–Volterra systems

titleCircuit design features of a stable two-cell system

*Supporting Information

§ INTRODUCTION

This Supporting Information is organized as follows. Section <ref> provides information about how we build LTI and LTV models of dynamics from time series data and previous perspectives on observability. Sections <ref> and <ref> provide mathematical formulations of Dynamic Sensor Selection and Structure Guided Sensor Selection respectively. Finally, Section <ref> outlines particular data sets and results obtained on the data. Supplementary figures are provided in Section <ref>, and the references in the Supporting Information are distinct from those in the main text.

§ RELATED WORK ON LEARNING DYNAMICS, OBSERVABILITY, AND SENSOR SELECTION

The problem of learning dynamics and selecting sensors or measurements is a wide field of study with a deep history <cit.>. Here, we clarify some of the particular algorithms we rely on for learning the dynamics of $\Avt{t}$ and survey classic measures of observability to distinguish our contributions. We consider the time-variant system with linear outputs

\begin{equation}\label{eq: SI model}
\begin{split}
\xvt{t+1} &=\Avt{t} \xvt{t}\textcolor{black}{+\Bv\uvt{t}}\\
\yvt{t} &=\Cvt{t} \xvt{t},
\end{split}
\end{equation}

where $\xvt{t}\in\R^n$ is the state vector, $\Avt{t}\in\R^{n\times n}$ is the state transition matrix, $\Cvt{t}\in\R^{p_t\times n}$ is the measurement matrix, and $\yvt{t}\in\R^{p(t)}$ is the system output. The measurements are lower dimensional than the state of our system (i.e. $p(t)\ll n$).

§.§ 1212Learning Dynamics

We focus on the parameterization of the state transition matrices $\Avt{t}$ from evenly spaced time series data $\Xv,$ where the data are of the form

\begin{equation*}
\Xv=\begin{bmatrix}
\xvt{0}&\xvt{1}&\dots&\xvt{T-1}&\xvt{T}\\
\end{bmatrix}.
\end{equation*}

Here, $\Xv\in\R^{n\times (T+1)}$ is a matrix with $n$ rows, where each row corresponds to a measurement, and $T$ columns, where each column contains all measurements at a particular instant in time. Consistent with many data, we assume the data are collected at even intervals in time (columns), and we will learn a model where each state variable comes from a single measurement (row). We first review Dynamic Mode Decomposition (DMD), then shift focus to the LTV systems under the Data Guided Control (DGC) model.

§.§.§ 1212Dynamic Mode Decomposition (DMD)

Here we outline the basic steps of DMD to provide a concise summary of the models we employ. For a comprehensive review on DMD and its applications in modeling biological and other complex systems, see <cit.> and references therein. Given time series data $\Xv,$ consider the first and last $T$ samples, denoted by $\Xv^-$ and $\Xv^+$ where

\begin{equation}
\begin{split}
\xvt{0}&\xvt{1}&\dots&\xvt{T-1}&\xvt{T}\\
\CodeAfter
\UnderBrace[yshift=25pt]{1-1}{2-4}{\Xv^-}
\OverBrace[yshift=3pt]{1-2}{1-5}{\Xv^+}
\end{NiceArray}
\text{ such that }
\Xv^-=\begin{bmatrix}
\xvt{0}&\xvt{1}&\dots&\xvt{T-1}\\
\end{bmatrix}\text{ and }\\    &\\
& \Xv^+=\begin{bmatrix}
\xvt{1}&\dots&\xvt{T-1}&\xvt{T}\\
\end{bmatrix}.
\end{split}
\end{equation}

The purpose of the DMD algorithm is to learn the best linear model of the dynamics $\xvt{t+1}=\Av\xvt{t}$ to explain the observed data. Learning the dynamics $\Av$ can be formalized by solving the matrix minimization

\begin{equation}\label{eq: DMD least squares problem}
\min_{A}\|\Xv^+-\Av \Xv^-\|^2_F,
\end{equation}

whose solution is given by $\Av=\Xv^{+}(\Xv^{-})^{\dagger}$, where $\dagger$ denotes the pseudo inverse. Consider the Singular Value Decomposition (SVD) $\Xv^{-}=\Uv\Sigma\Vv^{\top},$ where, $\Uv\in \R^{n\times r}$, $\Sigmam\in \R^{r\times r}$ and $\Vv\in \R^{(T-1) \times r}$ and $r\ll n$. The least squares problem <ref> has the solution

\begin{equation*}
\Av=\Xv^{+}(\Xm^{-})^{\dagger}=\Xm^{+}\Vv\Sigmam^{-1}\Uv^{\top},
\end{equation*}

which is the best linear model to explain the data $\Xv.$

Often, when we have multiple replicates of an experiment, rather than the data $\Xv$ being a $n\times T$ matrix, $\Xv$ will be a $n\times T\times r$ tensor, where $r$ denotes the number of replicates or separate instances of time series data generated from the system of interest. When there are multiple replicates, as in the case of all time series datasets we consider, we form $\Xv^+$ and $\Xv^-$ similarly by removing the first and last time points to generate $n\times (T-1)\times r$ tensors. Alternatively, the data can be averaged over multiple replicates to obtain a $n\times T$ matrix and be treated similarly as if only one replicate had been obtained.

§.§.§ 1212Model Reduction

DMD is classically applied to physical systems, which are often relatively low dimensional when compared with gene expression and biological systems.
The high dimensional number of variables in omics data lends itself to state space models of higher dimensions than most applications of DMD. For instance, the human genome contains on the order of $20$ thousand genes, each of which can be modeled as a state variable in $\xv$; however, such a model omits other important biological states or features such as isoform or protein expression. These large state space require model reduction to be computationally tractable, whereby the number of free parameters is reduced. Such an approach is consistent with many biological processes, where understanding the dynamics of only a few variables is sufficient to understand the dynamics of the system.

To study the dynamics of our high dimensional data $\Xv,$ we can change the coordinates to a lower dimensional space where it is easier to integrate, compute, and analyze our model.

Principal component (PC) based model reduction is standard to DMD. From computing the SVD of $\Xv^-$ in the DMD algorithm, the left singular vectors or PCs $\Uv$ map gene expression to principal component space.

To make such a change of coordinates, we apply the transformation $\xvttilde{t}=\Uv^\top\xvt{t}.$ Then the reduced state vector $\xvtilde$ will have $r$ components where $r$ is the rank of the data. Applying this transformation, we can produce a model reduced state transition matrix $\Amt,$ where $\Amt=\Uv^\top\Av\Uv.$ To see this reduction, consider

\begin{equation*}
\xvttilde{t+1}=\Um^{\top}\xvt{t+1}=\Um^{\top}\Av \xvt{t}=\Um^{\top}\Av \Um\xvttilde{t}=\Amt\xvttilde{t}.
\end{equation*}

From this, $\Amt$ can be derived directly from the data $\Xv$ as

\begin{equation*}
\Amt=\Um^{\top}\Av\Uv=\Uv^{\top}\Xv^{+}\Vm\Sigmam^{-1},
\end{equation*}

which is $r\times r$ matrix.

Returning to the observability of <ref>, we can cast the system output measurements $\yv$ in terms of the model reduced system as

\begin{equation*}
\begin{split}
% \nonumber to remove numbering (before each equation)
\xvttilde{t+1} &= \Amt\xvttilde{t}, \\
\yvt{t} &= \Cv \Uv \xvttilde{t},
\end{split}
\end{equation*}

where $\xvttilde{0}=\Uv^\top\xvt{0}.$
Thus, sensor selection problem can be cast in terms of reduced system.

§.§.§ 1212Switching Systems
Transitioning from LTI systems learned with DMD to LTV dynamics, consider a switching system, where there are two models of dynamics $\Av_1$ and $\Av_2$ that describe the flow of the system during two distinct phases in time. The dynamics of this model can be described as

\begin{eqnarray*}
\xvt{t+1} &=& \Av_{1} \xvt{t},\qquad 0\leq t \leq T_s\\
\xvt{t+1} &=& \Av_{2} \xvt{t},\qquad t> T_s,
\end{eqnarray*}

where $T_s$ denotes the time at which the dynamics switch from the first phase with $\Av_1$ dynamics to the second phase with $\Av_2$ dynamics. Dynamics such as these are well studied in social and communication systems and also characteristic of many biological dynamics such as bistable systems in cell regulation <cit.>.

Given data $\Xv_1$ is generated before $T_s$ and $\Xv_2$ after $T_s$, we can employ the DMD algorithm to model the dynamics $\Av_1$ and $\Av_2$ for each respective time period. Let $\Xv_i=\Uv_i\Sigmam_i\Vv_i^\top.$ Model reduction can then be applied to each time period, so that the dynamics of the system are expressed

\begin{eqnarray*}
% \nonumber to remove numbering (before each equation)
\xvttilde{t+1} &=& \Amt_1 \xvttilde{t},\qquad 0\leq t \leq T_s\\
\xvttilde{t+1} &=& \Amt_2 \xvttilde{t},\qquad t> T_s
\end{eqnarray*}

The full state space and reduced model are related to one another from the relation $\xvt{t}=\Uv_i^\top\xvt{t},$ where $i=1$ or $2,$ depending on the period of the system.

Given an initial condition $\xv_0$ the system evolves as

\begin{align*}
% \nonumber to remove numbering (before each equation)
\xvttilde{1} &= \Amt_1\Uv_1^{\top}\xvt{0}, &&\text{(model reduction)}\\
\xvttilde{2} &= \Amt_1\xvttilde{1}, &&\text{(step forward)}\\
\vdots&=\vdots&&\\
% \vdots &= \vdots \\
\xvttilde{T_s} &= \Amt_1\xvttilde{T_s-1}, &&\\%\text{(step forward)}\\
\xvttilde{T_s+1} &= \Amt_2\Uv_2^{\top}\Uv_1\xvttilde{T_s}, &&\text{(switch dynamics)}\\
\xvttilde{T_s+2} &= \Amt_2\xvttilde{T_s+1}, &&\text{(step forward)}\\
\vdots&=\vdots&&\\
% \vdots &= \vdots\\
\xvttilde{T} &= \Amt_2\xvttilde{T-1}. &&
\end{align*}

Thus, at any time point $t$, the reduced, one step, state transition matrix for a switching system can be written:

\begin{equation*}
\Amt_{t}=\begin{cases}
\Amt_1\Um_1^{\top}, & \mbox{if } t=0 \\
\Amt_1, & \mbox{if } 0<t\leq T_s \\
\Amt_2\Um_2^{\top}\Um_1, & \mbox{if } t=T_s+1 \\
\Amt_2, & \mbox{if } t>T_s+1 \\
\end{cases}.
\end{equation*}

§.§.§ 1212State Transition Matrices $\Phi$

Given time points time points $t_1$ and $t_2$ where $t_1\leq t_2$, we are interested in defining a single transition matrix between these times. For a LTI system, the transition matrix between $t_1$ and $t_2$ is found as $\Av^{t_2-t_1}$ or via a similar matrix exponentiation for the continuous case. For the LTV system <ref>, we can define a matrix $\Phi(t_2,t_1)$ that maps directly from time $t_1$ to $t_2$ as

\begin{equation*}
% \Phi(t_1,t_2)=\Amt_{t_2-1}\Amt_{t_2-2}\cdots \Amt_{t_1}.
\Phi(t_2,t_1)=\Avt{t_2-1}\Avt{t_2-2}\cdots \Avt{t_1}.
\end{equation*}

Then, given the initial conditions $\xv_{t_0}$ and the control signals $\uv_{t_0},\uv_{t_0+1},\cdots,\uv_{t-1}$ input to <ref>, the state $\xvt{t}$ can be written as

\begin{equation*}
\xvt{t}=\Phi(t,t_0)\xv_{t_0}+\sum_{i=t_0}^{t-1}\Phi(t,i+1)\Bvt{i}\uvt{i}.
\end{equation*}

The transition matrix $\Phi$ has following properties:

\begin{eqnarray*}
% \nonumber to remove numbering (before each equation)
\Phi(k_2,k_0) &=& \Phi(k_2,k_1)\Phi(k_1,k_0) \qquad k_0\leq k_1\leq k_2,\\
\Phi(k,k) &=& \In.
\end{eqnarray*}

These properties hold for the reduced transition matrix $\Tilde{\Phi},$ which is similarly defined as $\Tilde{\Phi}(t_1,t_2)=\Amtt{t_2-1}\Amtt{t_2-2}\cdots \Amtt{t_1}.$

§.§.§ 1212Data Guided Control (DGC) Model
The DGC model for approximating time variant linear systems was proposed by <cit.> to model the dynamics throughout cell reprogramming. This approach to learning dynamics considered the discrete-time LTV system with control:

\begin{equation}
\xvt{t+1}=\Avt{t}\xvt{t}+\Bv\uvt{t},
\end{equation}

where $\Bv_t$ is the control configuration matrix and $\uv$ is the input or control signal.

The fundamental assumption of the model is that gene expression of a population does not change considerably over time. Hence, the state transition matrix $\Av$ should be similar to the identity $\Iv.$ From this, the authors of <cit.> define $\Avt{t}$ as a rank one perturbation from the identity:

\begin{equation}
\Avt{t}=\Iv+\frac{(\xvt{t+1}-\xvt{t})\xvt{t}^\top}{\xvt{t}^\top\xvt{t}}.
\end{equation}

The authors then posed the challenge of selecting $\uvt{t}$ for cell reprogramming as an optimal control problem. Based on the success of this model for controlling cell dynamics, we chose to study its observability.

§.§ 1212Observability and Sensor Selection

The problem of observing dynamical systems is a broad area with deep roots in systems theory <cit.>. Mathematically, the ability to uniquely determine $\xv(0)$ is sufficient to call a system observable, since the knowledge of $\Av,$ $\xvt{0}$ and possibly a known control signal is sufficient to determine the state $\xvt{t}$ at any future point in time. Here we survey several classic tests for observability of linear systems, challenges associated with nonlinear observability, and how we can use observed outputs to estimate the full state of a system.

§.§.§ 1212Tests for Observability of Linear Systems

Classic tests for observability seek to provide binary criteria regarding if <ref> is observable. From the binary perspective, a system is mathematically observable if the full system initial condition $\xvt{0}$ can be uniquely determined.

In many instances, however, it is beneficial to consider more refined notions of observability, such as local observability, which is the ability to determine $\xvt{0}$ within a subset of possible states, or targeted observability, which is the ability to determine a subset of relevant state variables in $\xvt{0}$ <cit.>.

Here, we survey the Kalman Rank Condition, Popov-Belevitch-Hautus Test, and Structural Observability, three classical tests of observability.

1212Kalman Rank Condition.
The Kalman Rank condition for observability, likely the most famous criteria, guarantees a system is observable when the initial state of the system $\xv_0$ can be uniquely deteremined from the measurements $\yvt{0},\dots,\yvt{n-1}.$ As a test for this, the rank of the so called observability matrix $\mathcal{O}$ is compared with the dimension of the state $\xv$ so that when

\begin{eqnarray}
rank(\mathcal{O})=dim(\xv)\mbox{ where }\mathcal{O} = \begin{bmatrix}
\Cv\\
\Cv\Av\\
\Cv\Av^2\\
\vdots\\
\Cv\Av^{n-1}
\end{bmatrix}\mbox{ such that }    \begin{bmatrix}
\yvt{0}\\\yvt{1}\\\yvt{2}\\\vdots\\\yvt{n-1}
\end{bmatrix}=
\begin{bmatrix}
\Cv\\
\Cv\Av\\
\Cv\Av^2\\
\vdots\\
\Cv\Av^{n-1}
\end{bmatrix}\xvt{0}
\end{eqnarray}

has a unique solution for $\xvt{0},$ the system is observable. This approach is the foundation upon which subsequent tests and our estimation procedure are developed.

1212PBH Test.
Building upon the Kalman Condition, the Popov-Belevitch-Hautus (PBH) test, also known as the Hautus Lemma, guarantees observability where

\begin{eqnarray*}
rank(\mathcal{P})=\dim(\xv)\mbox{ where }
\mathcal{P}=\begin{bmatrix}
\mathbf{A} - \lambda \mathbf{I} \\
\mathbf{C}
\end{bmatrix},
\end{eqnarray*}

for all eigenvalues $\lambda$ of $\Av$ <cit.>. The rank deficiency of the PBH test can be seen as the existance of an eigenvector in the null space of the observability matrix $\mathcal{O}.$ By verifying $rank(\mathcal{P})=n$ for all $\lambda,$ the PBH test provides an equivalent check that $rank(\mathcal{O})=n$ or that $\xv(0)$ can be deteremined from the first $n$ measurements $\yvt{0},\dots,\yvt{n-1}$.

Although the Kalman Rank Condition and PBH Test offer convenient and consistent methods for assessing the observability of a LTI system, based on the matrices defining dynamics and observables, they are seldom utilized in practice due to three main challenges:

(C1) Identifiability: If we fail to identify $\Av$ to infinite precision, which will certainly be the case in any experimental system such as the ones studied here, then both the Kalman Rank Condition and PBH Test will almost certainly be observable <cit.>.

(C2) Cost: Explicitly constructing and determining the rank of $\mathcal{O}$ or $\mathcal{P}$ for all $\lambda$ is relatively expensive, particularly for high dimensional systems. For instance, many numerical schemes compute the SVD to determine the matrix rank, which is approximately a 3rd order polynomial operation, once the matrix is constructed.

(C3) Condition: Both the PBH and Kalman conditions provide a binary condition for observability where either YES a system is observable or NO it is not. In practice, however, many systems which are experimentally studied are unobservable, yet this does not deter us from gaining new insights. Rather, a graded degree of observability or targeted observability, where we are interested in only a few hidden states, is sufficient to understand internal dynamics from limited outputs.

To address C1-3, various alternatives have been proposed. For instance, Lin's structural controllability, in the next section, addressed both C1 and C2. Continuing this work, the measures of observability and formulations proposed in <ref> address C2 and C3 by allowing fast algorithms to provide scalar measures of observability, and the structural constraints of <ref> are designed to address C1.

1212Lin's Structural Observability.
To address the concern of system identifiability and, more recently, the cost involved in assessing observability for large or high dimensional systems, the theory of structural controllability was developed for LTI systems, and structural observability is formulated similarly <cit.>. To address (C1), Lin identified graph structures, based on the sparsity structures of $\Av$ and $\Cv$, that allow LTI dynamics to be observable based on the PBT test <cit.>. To address (C2), Liu, Slotine,  and Barabási developed the Minimum Inputs theorem from which a fast, graph based maximum matching algorithm can be applied to efficiently select sensor nodes.

Structural observability and controllability have been employed to study a variety of input/output systems across domains <cit.>. Yet there remain limitations to this approach as well <cit.>
For our work in particular, it leaves challenge C3 unanswered, as structural observability provides a binary YES/NO condition.

§.§.§ 1212Measures of Observability

The system output or observability energy provides a direct method to measure observability and address (C3). Output energy quantifies the amount of energy, defined as a norm, of the output measurements transmitted from a system $\yvt{0},\yvt{1},\dots$ with the equation

\begin{equation}\label{SI eq: energy}
\mathcal{E}=\sum_{t=0}\yvt{t}^\top\yvt{t}.
\end{equation}

When the output energy is small or zero, we lack useful information about the system's state. Therefore, one way to frame the sensor selection problem is as a maximization of the output energy.

The observability Gramian is a generalization of output energy from which many measures of observability have been proposed. Based on the relation that $\yvt{t}=\Cv\Av^t\xv(0),$ the Gramian generalizes <ref> as

\begin{equation}\label{eq: SI LTI G}
\G_o=\sum_{t=0}^\infty (\Av^t)^\top\Cv^\top\Cv\Av^t=\mathcal{O}^\top\mathcal{O}.
\end{equation}

This definition of the observability Gramian can also be derived from the Lyaponov equation $\Av^\top\G+\G\Av=-\Cv^\top\Cv.$ Due to the infinite summation, $\Av$ must be stable with all eigenvalues bound inside the unit circle to prevent $\G$ from diverging.

The eigenvalues of the $\G_o$ represent the relative observability of each system mode:

* The minimum eigenvalue $\lambda_{min}(\G_o)$ is a measure of the output energy for the least observable mode, and $\lambda^{-1}_{min}(\G_o)$ characterizes the  maximum estimation uncertainty.

* The maximum eigenvalue $\lambda_{max}(\G_o)$ is the measure of the output energy for the most observable mode, and $\lambda^{-1}_{max}(\G_o)$ characterizes the minimum estimation uncertainty.

* The eigenvector corresponding to the maximum eigenvalue, $\lambda_{max}(\G_o)$, is the direction with the largest gain (and thus most observable mode), therefore a small perturbation in that direction yields an output energy equivalent to that of a larger perturbation in the direction of the eigenvector corresponding to the minimum eigenvalue, $\lambda_{min}(\G_o)$ (direction of the least observable mode). Because the condition number $\kappa(\G_o) = \lambda_{max}(\G_o)/\lambda_{min}(\G_o)$ measures the ratio of maximum eigenvalue to minimum eigenvalue, an observability gramian $\G_o$ with a large condition number indicates that the output energy is dominated by some modes, while others are difficult to observe. Furthermore, $\kappa(\G_o^{-1})$ captures the shape of estimation uncertainty ellipsoid.

* The $\log det(\G_o)$ and equivalently $\log det(\G_o^{-1})$ measures the log of the volume of estimation uncertainty ellipsoid.

* The trace of the Gramian, i.e., $\textbf{tr}[\G_o]$ captures the average output energy, and equivalently $tr[\G_o^{-1}]$ measure the average estimation uncertainty.

Based on these, the following measures of observability have been proposed,

\begin{equation}\label{eq:measures}
\begin{split}
J_1(\G_o)=\textbf{tr}[\G_o^{-1}], \quad J_2(\G_o)=&\log det(\G_o^{-1}), \quad J_3(\G_o)=-\lambda_{min}(\G_o),\\ J_4(\G_o)=rank(\G_o), &\mbox{ and } J_5(\G_o)=\textbf{tr}[\G_o],
\end{split}
\end{equation}

which needs to be minimized to minimize the uncertainty in state estimate. Note that $J_1$ and $J_2$ are convex functions.  Also minimizing $\lambda^{-1}_{min}(\G_o)$ is equivalent to maximizing $\lambda_{max}(\G_o)$ which is a concave function. Hence, $J_3$ is also a convex function.

All the measures $J_i$ for $i=1,2,3$ are defined only when $\G_o$ is full rank. To handle such unobservable cases one could use metrics such as $\textbf{tr}[\G_o^{\dagger}]$, and corresponds to the average energy required to move the system around the observable subspace, or the log product of nonzero eigenvalues which relates to the “volume” of the subspace reachable with one unit of input energy.

§.§.§ 1212Least Squares Estimation

The mathematical notion of observability and our ability to estimate hidden states of a system from limited output are intimately related.

The relationship between these concepts begins with the Kalman Rank Condition, which guarantees mathematical observability for LTI systems

when there is a unique solution for $\xvt{0}$ to the system

\begin{equation}\label{SI EQ kalman 2}
\begin{bmatrix}
\yvt{0}\\\yvt{1}\\\yvt{2}\\\vdots\\\yvt{n-1}
\end{bmatrix}=
\begin{bmatrix}
\Cv\\
\Cv\Av\\
\Cv\Av^2\\
\vdots\\
\Cv\Av^{n-1}
\end{bmatrix}\xvt{0}.
\end{equation}

This problem is solved for $\xvt{0}$ directly as $\xvt{0}=\mathcal{O}^\dagger\Yv,$ where $\Yv$ is the left matrix in <ref> and $\mathcal{O}$ is the observability matrix as defined in <ref>.

This approach is readily adapted for LTV with control. Given a LTV system <ref>, with time dependent dynamics $\Avt{t}$, measurements $\Cvt{t}$, input matrices $\Bvt{t}$, and control signals $\uvt{t}$, the output measurements of the system evolve over time as

\begin{eqnarray*}
\yvt{0} &=& \Cvt{0} \xvt{0}, \notag \\
\yvt{1} &=& \Cvt{1} \xvt{1}=\Cvt{1}\Phi(1,0)\xvt{0}+\Cvt{1} \Bvt{0}\uvt{0} ,\notag\\
\vdots &=& \vdots \notag\\
\yvt{t} &=& \Cvt{T} \xvt{t}=\Cvt{T}\Phi(T,0)\xvt{0}+\Cvt{T} \sum_{i=0}^{T-1}\Phi(T,i+1)\Bvt{i}\uvt{i}. \notag
\end{eqnarray*}

Thus, one can express the system as,

\begin{equation}\label{eq:ls}
\Yv=\mathcal{O} \xvt{0}+\Uv\qquad\text{ with the solution }\qquad\xvt{0}=\mathcal{O}^\dagger(\Yv-\Uv)
\end{equation}

\begin{equation}\label{eq:}
\Yv=\begin{bmatrix}%{c}
\yvt{0} \\
\yvt{1}\\
\vdots \\
\yvt{T}
\end{bmatrix},%\right)
\quad\mathcal{O}=\begin{bmatrix}%{c}
\Cvt{0}  \\
\Cvt{1}\Phi(1,0)\\
\vdots \\
\Cvt{T}\Phi(T,0)
\end{bmatrix},\quad \text{and} \quad\Uv=
\begin{bmatrix}
\mathbf{0} \\
\Cvt{1} \Bvt{0}\uvt{0}\\
\vdots \\
\Cvt{T} \sum_{i=0}^{T-1}\Phi(T,i+1)\Bvt{i}\uvt{i}\\
\end{bmatrix}.
\end{equation}

There is a unique solution for $\xvt{0}$ in the above two equations when $\mathcal{O}^\dagger$ is full rank, motivating the Kalman Rank condition to check that $rank(\mathcal{O})=dim(\xv).$

Assuming that the measurements $\yv$ are corrupted by independent and identically distributed (i.i.d.) zero-mean Gaussian noise with normal distribution $\mathcal{N}(0,\sigma \In)$, using a weighted least squares formulation, the minimum variance estimate for $\hat{\xv}^*(0)$ is given by,

\begin{equation}\label{eq: LS estimation of xhat}
\hat{\xv}^*(0)=(\mathcal{O}^\top\mathcal{O})^{-1}\mathcal{O}^\top (\Yv-\Uv).
\end{equation}

The estimate of $\hat{\xv}_0^*$ has error covariance $\Pv_{\xv_0}$

\begin{equation*}
\Pm_{\xvt{0}}=\mathbb{E}\big[(\hat{\xv}^*(0)-\xvt{0})(\hat{\xv}^*(0)-\xvt{0})^\top\big]=\sigma(\mathcal{O}^\top\mathcal{O})^{-1}=\sigma (\G_o)^{-1}.
\end{equation*}

The least squares estimator is known to be efficient (i.e., the Cramer-Rao lower bound is achieved), therefore the Fisher Information Matrix (FIM) $\Fm =\sigma^{-1}\mathcal{O}^\top\mathcal{O}$ is exactly the inverse of the estimation covariance <cit.>.

To see the relation between the LS estimation proceedure and the Gramian, note that the observability Gramian can be expressed in terms of $\mathcal{O}$ based on the relationship

\begin{eqnarray}
\mathcal{O}^\top\mathcal{O} &=&
\begin{bmatrix}
\Cvt{0}^\top & \Phi(1,0)^\top\Cvt{1}^\top & \cdots & \Phi(T,0)^\top\Cvt{T}^\top \\
\end{bmatrix}
\begin{bmatrix}
\Cm_0  \\
\Cm_1\Phi(1,0)\\
\vdots \\
\Cm_T\Phi(T,0)
\end{bmatrix},\notag\\
&=&\sum_{i=0}^{T}\Phi(i,0)^\top\Cvt{i}^\top\Cvt{i}\Phi(i,0)=\G_o(0) \label{eq:invO}.
\end{eqnarray}

By relation (<ref>), the observability gramian  $\G_o$ is proportional to the FIM and inversely proportional to the estimate covariance. Consequently,  the eigenvalues of the observability Gramian directly control the Fisher information and inversely control the estimation covariance, and can be used to define sensor selection measures/metrics. Moreover, when $\G_o$ is full rank, the least squares estimation of $\hat{\xv}(0)^*$ in <ref> is well determined.

Despite this relationship between observability and estimating $\xvt{0},$ to date, the marriage of these concepts for biological systems has been prevented due to practical considerations.

* Noise: Experimental observations or measurements in a real system are corrupted with noise.

* Scalability: the least squares approach requires a significant number of observations to achieve the theoretical guarantee of observability.

In the experimental setting, our measurements $\yvt{0},\yvt{1},\dots$ will always be corrupted with noise, and for high dimensional systems where $\xv\in\R^n$ and $n$ is large, such as genomics data ($n\gg 20,000$), it is cost prohibitive and experimentally infeasible to obtain $n$ measurements. Fortunately, through careful sensor selection truncated observability matrices $\mathcal{O}$, where the number of rows are truncated to only include time points for which we have data, can be designed so that least squares estimation is still well conditioned. We aim to address these considerations in addition to C1-3 in the following sections regarding Dynamic and Structure Guided Sensor Selection.

§ DYNAMIC SENSOR SELECTION

In this section, we return to <ref> and address the challenge of selecting the best measurement matrices $\Cvt{t}$ when given known dynamics and state transition matrices $\Avt{t}.$ We first review some basic objects and control perspectives on observability for LTI systems and then present two formulations for sensor selection for time variant systems.

§.§ 1212Oscillator Networks

Oscillator networks are characterized by internal dynamics of each oscillator, as well as interactions among them. This is described by the equation:

\begin{equation}\label{SI eq oscillator}
\dfrac{d\xvt{t}}{dt} = \textbf{F}(\xvt{t}, \mu) - \textbf{L}(\xvt{t}),
\end{equation}

as seen in <ref>.
Here, $\xv$ is a vector representing the locations or values of each oscillator, $F$ is the dynamics of individual oscillators with internal parameters $\mu,$ and $L$ is the diffusion operator specific to the network structure.

The dynamics of $\Fv(\xvt{t})$ can be any oscillator, such as the Van der Pol or Andorov-Hopf oscillators, where each oscillator may possess its unique parameters. Additionally, the coupling between oscillators in $\mathbf{L}(\xvt{t})$ can occur across networks with diverse directionalities and weights. This framework via <ref> encompasses several well-known oscillators, including: (1) the Kuromoto oscillator, in which individual oscillators positioned along the unit circle interact with sinusoidal functions <cit.>; (2) Turing's Equation, which generates oscillations by coupling of two stable systems, and Smale's related work <cit.>.

Continuing the work of Kuromoto, Turing, and Smale, systems of this form have garnered recent interest in the study of higher order structures. For instance, network motifs couple groups of oscillators to form complex and emergent behavior that have been studying in network science <cit.>, systems biology <cit.>, communication systems <cit.>, among other systems. In network biology, genetic circuits such the Repressilator <cit.>, Goodwin Oscillator <cit.>, Toggle Switch <cit.> can be mathematically expressed in the form of <ref>.

§.§.§ 1212Repressilator

at (0,0) ;
at (-5.5,-2) 1616$\begin{aligned}
\dfrac{dm_i}{dt} & = -m_i+\dfrac{\alpha}{1+p_j^n}+\alpha_0\\
\dfrac{dp}{dt} & = -\beta(p_i-m_i)
\end{aligned}$

Repressilator. (A) The Repressilator contains three nodes, each of which represents the expression of a gene. These nodes are coupled in an inhibitory chain where as the state or concentration of one gene increases, it will inhibit expression of the next gene in the chain. (B) The dynamics of a single node are expressed with these differential equations where $m_i$ and $p_i$ denote the concentrations of mRNA and protein associated with gene $i$ respectively. Consistent with the biology, production of protein is positively regulated by the mRNA concentration and both mRNA and protein negatively regulate themselves. Additionally, the production of mRNA of gene $i$ decreases with the concentration of protein for gene $j$ to encode the inhibitory relationship. The full system of equations for all three genes together is in <ref>(C) We simulate several oscillations of the Repressilator from random initial conditions with $\alpha=20,\alpha_0=0,\beta=2$ and $n=2.1.$

The canonical Repressilator shown in <ref> is one such instance of a network of coupled oscillators, each of which represents the expression from DNA to mRNA and finally protein. In Escherichia coli., the expression from mRNA to protein of lacI, tetR, and cl can be seen as three osicillators, where the regulation the mRNA and protein concentration for a single gene is both regulated internally ($\Fv(\xvt{t})$) and repressed by coupled proteins ($\Lv(\xvt{t})$). Borrowed from <cit.>, the continuous deterministic Repressilator has dynamics

\begin{eqnarray*}
% \overbrace{\dfrac{dm_{lacI}}{dt}}^{\dfrac{d\xv}{dt}}
\dfrac{dm_{lacI}}{dt} &= \overbrace{-m_{lacI}+\alpha_0\phantom{\dfrac{x}{x}}}^{\Fv(\xv)}&\overbrace{+\dfrac{\alpha}{1+p_{cl}^n}}^{-\Lv(\xv)},\\
\dfrac{dm_{tetR}}{dt} &= -m_{tetR}+\alpha_0&+\dfrac{\alpha}{1+p_{lacI}^n},\\
\dfrac{dm_{cl}}{dt} &= -m_{cl}+\alpha_0&+\dfrac{\alpha}{1+p_{tetR}^n},\\
\end{eqnarray*}

\begin{equation}\label{SI: rep full guy}
\begin{split}
% \overbrace{\dfrac{dp_{lacI}}{dt}}^{\dfrac{d\xv}{dt}}
\dfrac{dp_{lacI}}{dt} &= \overbrace{-\beta p_{lacI}\phantom{\dfrac{x}{x}}}^{\Fv(\xv)}\overbrace{+\beta m_{lacI}\phantom{\dfrac{x}{x}}}^{-\Lv(\xv)}, \\
\dfrac{dp_{tetR}}{dt} &= -\beta p_{tetR}+\beta m_{tetR}, \\
\dfrac{dp_{cl}}{dt} &= -\beta p_{cl}+\beta m_{cl}.
\end{split}
\end{equation}

Here, there are have six chemical species describing the mRNA ($m$) and protein ($p$) concentrations of the lacI, tetR, and cl genes; $\alpha, \beta,$ and $n,$ are parameters; and there is a pair of equations to describe the dynamics of mRNA and protein concentrations separately. Each form of the equations can be expressed in terms of <ref>, where $\Fv(\xvt{t})$ describes how the concentrations change as a function of the concentration of a given species, and $\Lv(\xvt{t})$ describes how the concentration changes as a function of either the upstream mRNA or the repressing protein.

§.§.§ 1212Oscillator Types

Oscillator networks, particularly those found in gene regulatory or neurological systems, can be directly modeled from data, as we do in the analysis of various genomics and EEG time series data; however, in addition to empirical approaches, we may also explore theoretical models of such oscillators. We turn to the following well known oscillators:

* Van der Pol (VP):

\begin{equation}\label{eq: VP and AH oscillators}
\dfrac{d^2x}{dt^2} - \mu (1 - x^2) \dfrac{dx}{dt} + x = 0 \qquad\text{ or }\qquad
\begin{cases}
\dfrac{{dx}}{{dt}} &= \mu \cdot (1 - y^2) \cdot x - y\\
\dfrac{{dy}}{{dt}} &= x
\end{cases}
\end{equation}

* Goodwin (GW):

\begin{equation}\label{eq: goodwin}
\begin{cases}
\dfrac{dx_1}{dt} & = \dfrac{\alpha}{1+x_3^n}-x_1\\
\dfrac{dx_2}{dt} & = x_1-x_2\\
\dfrac{dx_3}{dt} & = x_2-x_3\\
\end{cases}
\end{equation}

* Andronov-Hopf (AH):

\begin{equation}\label{eq: Andronov-hopf}
\begin{cases}
\dfrac{dx}{dt}&= ax-by-x(x^2+y^2)\\
\dfrac{dy}{dt}&= bx+ay-y(x^2+y^2)
\end{cases}
\end{equation}

The VP oscillator was introduced to model vacuum tubes by the Dutch physicist Balthasar van der Pol (1889-1959).

This second order differential equation and its representation as a pair of first-order differential equations via the Koopman theory, is a non-conservative oscillator that exhibits nonlinear damping.

This model has since been applied to study several other biological and complex systems, such as neurological circuit action potentials, and has been well studied <cit.>.

The GW oscillator, proposed by biologist Brian Goodwin (1931-2009), was an early model of a genetic oscillator with three variables to model concentrations of RNA, proteins, and a final product all produced from the same gene and regulated with negative feedback <cit.>. During the early to mid-1960s, this model of a genetic oscillator was proposed shortly after François Jacob and Jacques Monod introduced their model of gene regulation <cit.>. The oscillator was developed as a Hamiltonian system and an early proponent of using Hill functions in biology.

The AH oscillator, after the physicist Aleksandr Andronov (1901-1952) and astronomer and mathematician Eberhard Hopf (1902-1983), as shown in <ref> is one particular form of a general AH oscillator, which is a pair of coupled oscillators $x$ and $y$ that exhibit a Hopf bifurcation (when a critical point becomes a limit cycle); the Hopf bifurcation is sometimes also referred to as the Poincaré-Andronov-Hopf bifurcation to provide attribution to Poincaré and Andronov for their work studying this bifurcation <cit.>. Numerous other ecological and biological oscillators, such as the Lotka–Voltera, SIR, Hodgkin–Huxley models, exhibit Hopf bifurcations as well <cit.>.

1212Internal Oscillator Parameters.
The trajectories an oscillator networks are governed in part by the internal dynamics of each oscillator and governed by parameters associated with individual oscillators.

The VP oscillator is governed by the parameter $\mu$, which controls the stability and limit cycle of the system. In <ref>.A, the dynamics of three instances of the VP oscillator are shown. As $\mu$ changes from less than zero to greater than zero, the limiting behavior of the system changes from the origin being a stable critical point to a limit cycle with fast and slow transitions. The top row in <ref>.A plots the position $x$ in terms of time, and the stability of this system is interpreted based upon how the oscillation amplitud changes between periods. For instance, in the top row where $\mu<0$, the contraction of the position over time indicates stability. The bottom row of in <ref>.A plots the position $x$ relative to the velocity $y$, and there the stability can be seen based upon the shift of the trajectory toward the origin or in a stable limit cycle. For instance, in the bottom row where $\mu>0,$ the trajectory enters a limit cycle where it remains throughout time, similar to the fixed amplitude oscillations above. By varying $\mu,$ the behavior of individual oscillators is changed and can effect the dynamics of an oscillator network.

The GW oscillator of <ref> is governed by two parameters $\alphav$ and $n.$ Often, and in the original parlance of Goodwin, this oscillator is written with more parameters, but the number of free parameters may be reduced to 2 in the case of equal degredation rates when written in the dimensionless form, which is useful for the purpose of considering its birfurcation behavior <cit.>. Regardless, the Goodwin oscillator exhibits a Hopf bifurcation, as illustrated in <ref>.B. By fixing the value of $n,$ and varying $\alpha,$ we see that for certain values of $\alpha$ the system is stable while other values produce internal oscillations in the form of a limit cycle. The AH, VP, Liénard, and many other oscillators exhibit Hopf bifurcations as well.

(A) The dynamics of a single VP oscillator depend on its parameter $\mu.$ The top row shows the time parameteric state of the oscillator where $x$ from the left of <ref> is shown as a function of time and the bottom row plots the position $x$ relative to the velocity $y$ based on the right side of <ref>. (B) We fix $n=20$ and vary the bifurcation parameter $\alpha$ from 0.5 to 2. For each value of $\alpha,$ we simulate a random trajectory, shown in grey, of the oscillator, and plot the end point of the trajectory in black. On the left hand side of this plot, all the trajectories arrive at a final point near each other indicating stability as a fixed point, and on the right hand side, the black dots and end of the trajectories are not necessarily near each other but rather positioned along a limit cycle. This transition occurs near 1.6 and demonstrates the existance of a Hopf bifurcation. Similar examples of the Hopf bifurcation on the AH, VP, and Liénard system are shown below.

§.§.§ 1212Diffusion and Strong Coupled Oscillators

The trajectories of oscillator networks are governed in part by the network structure and external interactions between oscillators.

Diffusive coupling of oscillators forces the state of an oscillator toward the state of its neighboring osillators. In particular, when the oscillators are coupling on a graph $\mathcal{G}$ with a Laplacian $\Lv,$ the diffusive coupling of a single oscillator $\mathbf{x}_i$ can be written as

\begin{equation}
\Lv(\mathbf{x}_i) = \lambda\sum_{j\in\mathcal{N}(i)}\mathbf{x}_i-\mathbf{x}_j,
\end{equation}

where $\mathcal{N}(i)$ is the neighbors or oscillators that are adjacent to oscillator $i$ on the graph $\mathcal{G},$ $\lambda$ denotes the coupling strength, and the vector $\mathbf{x}_i$ represents all state variables of the oscillator.

Returning to <ref>, there are three VP oscillators coupled in a directed chain. Each of these oscillators has the internal dynamics $\Fv(\xv_i)$ of <ref>, and the chain of interactions in <ref>.D defines the graph laplacian $\Lv(\xv)$ two edges, from oscillators $\xv_1$ to $\xv_2$ and from $\xv_2$ to $\xv_3.$ The full dynamics of this system are written as

\begin{eqnarray}
\dfrac{{dx_1}}{{dt}} =& \overbrace{\mu_1 \cdot (1 - y_1^2) \cdot x_1 - y_1}^{\Fv(\xv,\mu)}&\overbrace{-0\phantom{(^2}\ \quad\quad\quad,}^{-\Lv(\xv)}\notag\\
\dfrac{{dy_1}}{{dt}} =& x_1&-0\quad\quad\quad\quad,\notag\\
\dfrac{{dx_2}}{{dt}} =& \mu_2 \cdot (1 - y_2^2) \cdot x_2 - y_2&-\lambda(x_2-x_1),\label{eq full dynamics of figure 1}\\
\dfrac{{dy_2}}{{dt}} =& x_2&-\lambda(y_2-y_1),\notag\\
\dfrac{{dx_3}}{{dt}} =& \mu_3 \cdot (1 - y_3^2) \cdot x_3 - y_3&-\lambda(x_3-x_2),\notag\\
\dfrac{{dy_3}}{{dt}} =& x_3&-\lambda(y_3-y_2).\notag
\end{eqnarray}

Here, there are three Van der Pol oscillators written in their two dimensional form $\Fv(\xv,\mu)$, each with its on internal state $\xv_i=\begin{bmatrix}
\end{bmatrix}^\top$, and internal parameters $\mu_i$ discussed in the prior section. In eq. <ref>, $x_1$ and $y_1$ are not influenced by external parameters from other oscillators, so the long term behavior of the oscillator network is controlled by oscillator 1.

From Lin's perspective of structural observability/control, oscillator 1 is called a root node. This means all other nodes can be influenced from oscillator one and it is uneffected by any node. As a result, from the perspective of structural observability/control, oscillator 1 is a good sensor/actuator. This is based on the notion that the LTI system defined on the adjacency matrix of the graph from <ref>.D is observable/controllable with oscillator 1 as a sensor/actuator according to the PBH test.

One factor not considered by Lin's structural perspective is the strength or weight of the coupling between oscillators. In eq. <ref>, all oscillators are coupled with an identical parameter $\lambda>0,$ but in practice, each interaction could have its own coupling strength. Since coupled oscillators interact by adjusting their phase relative to one another, the coupling coefficients $\lambda,$ control both the speed and direction with which oscillators interact <cit.>.

In <ref> and <ref>, for instance, several instances of coupled VP and AH oscillators are shown. We consider oscillator networks containing either 2 or 5 oscillators coupled in a ring or cycle network where all interactions have a weight $\lambda.$ In the first two rows, we consider the case where $\lambda_1=1$ and $\lambda_2=-1$ such that if the oscillators were uncoupled, one would decay quickly toward the origin and one would exhibit a limit cycle. In the top row, $\lambda>0$ such that the oscillator in the limit cycle induces synchronized oscillations in the otherwise stable oscillator, and in the second row, where $\lambda<0,$ so that the induced oscillations of the stable VP oscillator are out of phase with the oscillator in the limit cycle. Rows 3 and 4 recreate rows 1 and 2 where both oscillators are parameterized with $\mu>0$ so that both oscillators exhibit limit cycles. Here the coupling of $\lambda$ synchronizes the oscillations rather than induces new ones, and comparing between the left and right columns where the magnitudes of $\lambda$ change, we see that as the coupling strength increases, the oscillators shift into synchrony or phase locking quicker. Finally, rows 5 and 6 recreate identical results as in rows 3 and 4, with an expanded number of oscillators so that the synchronization and phase locking is more clear.

Coupled VP oscillators are used in the above 12 simulations. The oscillators are coupled in rings, and rows 1-4 have networks of 2 oscillators whereas rows 5 and 6 have 5 oscillators. Each row illustrates one type of dynamics, where the left column illustrates weak coupling between the oscillators and the right column illustrates strong coupling, governed by the parameter $\lambda.$ Row 1: an oscillator with $\mu=1$, which alone exhibits a limit cycle, is coupled with an oscillator where $\mu=-1$ that would decay. The oscillating node induces synchronized oscillations in the network. Row 2: Similar set up as row 1, now with $\lambda>0$ so that phase locking is induced. Row 3: synchronization between 2 oscillators nodes. Row 4: Phase locking between 2 oscillating nodes. Row 5: synchronization between several oscillators. Row 6: phase locking between several oscillators.

Coupled AH oscillators are used in the above 12 simulations. Whereas VP oscillators have a single position variable $x,$ AH oscillators have two variables $x$ and $y.$ The behavior is qualitatively similar between the $x$ and $y$ variables, so only the position of the $x$ variables for each oscillator is shown in these simulations. The interpretation of these twelve plots mirrors that of the corresponding simulations of VP oscillators detailed in <ref>

§.§ 1212Sensor Selection Formulations

Based on the oscillator networks as a guiding example and the reviewed methods of observability in <ref>, here we propose two methods for dynamically selecting sensors.

Consider the system

\begin{equation}\label{eq: SI LTI}
\begin{split}
\begin{cases}
\xvt{t+1}&=\Avt{t}\xvt{t}\\
\yvt{t}&=\Cvt{t}\xvt{t},
\end{cases}
\end{split}
\end{equation}

which is equivalent to <ref> without control. We summarize two methods of sensor selection and designing $\Cv$ for <ref>, which are extensions of similar methods for LTI systems.

* Energy based selection: One graded approach to maximize observability is to maximize the output energy of the system, formalized as:

\begin{equation}\label{SI eq energy problem}
\max_{\Cvt{t}}\mathcal{E}=\max_{\Cvt{t}}\xvt{0}^\top\G\xvt{0}\text{ subject to }\Cvt{t}\Cvt{t}^\top=\Iv.
\end{equation}

A similar formulation with LTI dynamics with fixed outputs was used by Hasnain et. al. to select biomarkers <cit.>. This maximization problem is a quadratic program with linear constraints, and particularly when $\Tilde{\Av}$ has been computed from DMD, the sensor rankings from this optimization can be obtained quickly.

* Gramian based selection: Sensors can be selected to optimize a particular function of the Gramian that corresponds to the observability of the system. We consider the maximization

\begin{equation}\label{SI eq gramian problem}
\max_{\Cvt{t}} J(\G),
\end{equation}

where $J$ is one of the five functions discussed in <ref>.
Summers et. al. utilized this formulation for fixed sensor selection on LTI systems <cit.>. The five proposed functions $J$ of the Gramian are submodular, making them conducive to greedy algorithms for sensor selection, and when the trace is used, <ref> is solved with a linear program.

Since the objective functions of the optimization problems in <ref> and <ref> are continuous, (C3) is addressed by providing scalar value of observability rather than a binary criteria.

Moreover, since both optimizations allow for the use of model reduction, and can be solved with fast algorithms to address (C2).

In the subsequent sections, we provide further detail regarding how to solve these optimization problems and select sensors.

§.§ 1212Output Energy Maximization

Here we provide our method to solve <ref> based on its Lagrangian dual form. We first discuss how this problem is solved when the sensors $\Cv$ are fixed for all time and then consider the dynamic selection of sensors.

§.§.§ 1212Time Invariant Sensors

The objective is to select sensors $\Cv$ that maximize the signal or output energy of the system $\mathcal{E}$ where

\begin{equation*}
\mathcal{E}=\sum_{i=t_0}\yvt{i}^\top\yvt{i}=\sum_{i}\xvt{0}^\top\Phi(i,t_0)^{\top}\Cv^\top\Cv\Phi(i,t_0)\xvt{0},
\end{equation*}

which is formalized as

\begin{equation}\label{SI eq cv LTI}
\max_{\Cv}\mathcal{E}\text{ subject to }\Cv^\top\Cv=\Iv.
\end{equation}

<Ref> is the fixed sensor formulation of <ref>. Hasnain et. al. shows that
the Lagrangian dual formulation of this problem is

\begin{equation*}
\max_{\Cv}\mathcal{E+L}\text{ where }\mathcal{L}=\textbf{tr}((\Cv\Cv^\top-\Iv)\Dv),
\end{equation*}

where $\Dv$ are the dual variables <cit.>. Following eq. 5 of <cit.>,

\begin{equation}
\dfrac{\partial(\mathcal{E+L})}{\partial\Cv^\top}=2\Gv\Cv^\top-2\Cv^\top\Dv=0,\text{ such that }\Gv\Cv^\top=\Cv^\top\Dv.
\end{equation}

This last expression implies that the eigenvectors of $\Gv$ are the sensor weights or the importance of each sensor at a critical point of the signal output energy with respect to the sensors.

We extend the approach to the selection of time varying sensors.

§.§.§ 1212Time Variant Sensors

Suppose we seek to maximize the signal energy at time $t$ such that our objective function is $\mathcal{J}_t$ where

\begin{equation}
\mathcal{E}(t)=\yvt{t}^\top\yvt{t}=\xvt{0}^\top\Phi(t,t_0)^\top\Cvt{t}^\top\Cvt{t}\Phi(t,t_0)\xvt{0}.
\end{equation}

Applying the prior Lagrangian formulation here, we seek to maximize $\mathcal{E}(t)$ with respect to $\Cvt{t}$ at all times. Since $\Cvt{t_{1}}$ and $\Cvt{t_2}$ are independent of one another for all $t_1$ and $t_2,$ we can maximize the energy at each time $\mathcal{E}(t)$ independently of one another. To do so, we extend the fixed sensor selection optimization of $\mathcal{E}$ to an optimization of $\mathcal{E}(t),$ using a similar approach as in <cit.>. This is formulated as follows:

\begin{equation}\label{SI eq lagrangian dual form for energy}
\begin{split}
\dfrac{\partial(\mathcal{E}_t+\mathcal{L})}{\partial\Cvt{t}^\top} & =\dfrac{\partial}{\partial\Cvt{t}^\top}\bigg(\xvt{0}^\top\Phi(t,t_0)^\top\Cvt{t}^\top\Cvt{t}\Phi(t,t_0)\xvt{0}-\textbf{tr}((\Cvt{t}\Cvt{t}^\top-\Iv)\Dv)\bigg)\\
& =\dfrac{\partial}{\partial\Cvt{t}^\top}\bigg(\textbf{tr}\big(\xvt{0}^\top\Phi(t,t_0)^\top\Cvt{t}^\top\Cvt{t}\Phi(t,t_0)\xvt{0}\big)-\textbf{tr}((\Cvt{t}\Cvt{t}^\top-\Iv)\Dv)\bigg)\\
& =\dfrac{\partial}{\partial\Cvt{t}^\top}\bigg(\textbf{tr}\big(\Cvt{t}\Phi(t,t_0)\xvt{0}\xvt{0}^\top\Phi(t,t_0)^\top\Cvt{t}^\top\big)-\textbf{tr}((\Cvt{t}\Cvt{t}^\top-\Iv)\Dv)\bigg).\\
\end{split}
\end{equation}

At this point, let

\begin{equation*}
\Gv(t, t_0)= \Phi(t,t_0)\xvt{0}\xvt{0}^\top\Phi(t,t_0)^\top.
\end{equation*}

Then, returning to <ref> we can obtain the optimal sensors $\Cvt{t}$

\begin{equation}
\begin{split}
\dfrac{\partial(\mathcal{E}_t+\mathcal{L})}{\partial\Cvt{t}^\top} & =\dfrac{\partial}{\partial\Cvt{t}^\top}\bigg(\textbf{tr}\big(\Cvt{t}\Gv(t,t_0)\Cvt{t}^\top\big)-\textbf{tr}((\Cvt{t}\Cvt{t}^\top-\Iv)\Dv)\bigg)\\
& = 2\Gv(t,t_0)\Cvt{t}^\top-2\Cvt{t}^\top\Dv=0.
\end{split}
\end{equation}

This has a similar interpretation to the time invariant case, were the eigenvalues of $\Gv(t,t_0)$ denote the contribution of each state variable to observability at time $t.$

To solve this and select sensors from time $t_0,\dots,t,$ we must form $\Gv(t, t_0)$ for all $t.$ This requires integrating the system forward from the initial conditions $\xvt{0},$ which can be performed efficiently using model reduction, and then computing the largest eigenvector of the matrices $\Gv(t, t_0)$, for which there are fast algorithms.

§.§ 1212Gramian Based Observability

Here we provide our method to solve <ref> with standard optimization techniques. We begin by highlighting a few key properties of the observability Gramian for LTV systems, then provide the integer programming formulation of <ref> as an integer programming. Finally, we provide a continuous relaxation of the integer programming problem, making it solvable as a linear program.

§.§.§ 1212Integer Programming Formulation

To formulate the sensor selection problem in terms of integer programming, consider the case where there is a binary variable denoting whether or not each state variable $x_i\in \xv$ is observed or measured at time $t.$ Let $\cv_{tj}$ be the $j$-th row of the matrix $\Cm_t \in R^{n\times p_t}$. Then the LTI observability Gramian can be reweitten as

\begin{eqnarray*}
% \nonumber to remove numbering (before each equation)
\G_o &=& \sum_{i=0}^{}\Phi(i,0)^\top\Cvt{i}^\top\Cvt{i}\Phi(i,0), \\
&=&  \sum_{i=0}^{} \sum_{j=1}^{p_k}\Phi(i,0)^\top(\cv_{ij})^\top\cv_{ij}\Phi(i,0).
\end{eqnarray*}

Let, $\alpha_{tj}\in\{0,1\}$ be a binary variable, which indicates whether $x_j$ variable is measured at time $t$. Then one can express above relation, as

\begin{eqnarray*}
% \nonumber to remove numbering (before each equation)
\G_o(\alphav) &=&\sum_{i=0}^{} \sum_{j=1}^{n}\alpha_{ij} \Wm_{ij},
\end{eqnarray*}

where, $\alphav=(\alpha_{11},\alpha_{21},\cdots, \alpha_{(T+1)n})^\top\in R^{n(T+1)},$ where $T$ is the upper bound on the first summation used to compute the observability Grammian, and

\begin{equation*}
\Wm_{ij}=\Phi(i,0)^\top(\cv_{ij})^\top\cv_{ij}\Phi(i,0).
\end{equation*}

Within this representation, note that $\cv_{ij}$ is a row vector with $j$-th entry as $1$ and zero otherwise. Note that each $\Wm_{ij}$ is a positive semidefinite matrix. Because composition with affine mappings preserves convexity, each measure $J_i(\G_o)$ discussed in the Section <ref> is also a convex function $J_i(\G_o(\alphav))$ of the sensor selection variable $\alphav$.

Given the sensor selection variables $\alphav$, the sensor selection problem can be written as a mixed-integer convex problem,

\begin{equation}\label{eq:optint}
\min_{\alphav}  J(\G_o(\alphav)) \quad \mbox{subject to} \quad \sum_{j=1}^{n}\alpha_{tj}\leq p_t,\quad
t=0,1,\cdots \quad \text{where }\quad 0\leq \alphav_{kj}\in\{0,1\}.
% && \alphav_{kj}\in\{0,1\}.  \notag
\end{equation}

Here, the first constraint restricts number of selected sensors to be no more than $p_t$ for each time point. Because mixed-integer programs do not scale well for large problems, a convex relaxation to (<ref>) provides a useful solution alternative.

§.§.§ 1212Continuous Relaxation

In the continuous relaxation, the observation of variable $j$ at time $t$ is relaxed to the interval $\alpha_{ij}\in[0,1].$ This leads to the convex program,

\begin{equation}
\begin{split} \label{eq:optrelax}
\min_{\alphav}  J(\G_o(\alphav)) \quad
\mbox{subject to}\quad\sum_{j=1}^{n}\alpha_{tj}\leq p_t,\quad t=0,1,\cdots,\quad \text{where }\quad 0\leq \alphav_{tj}\leq 1.
\end{split}
\end{equation}

The advantage of the relaxation is that it can be solved in time that is polynomial in the number of variables using efficient techniques such as interior point methods. Furthermore, if the solution to the relaxed problem is such that $\alpha_{ij}\in \{0,1\}$ (within numerical tolerance), then the original mixed-integer problem has been solved. The relaxation serves two roles — an approximate (suboptimal) solution to the mixed-integer problem by rounding $\alpha_{ij}$, and, in some cases, a fast optimal solution to the mixed-integer problem.

In both the mixed-integer problem (<ref>) and the convex relaxation (<ref>), the desired number of sensors was explicitly constrained to be $p_t$. Another approach is to allow the number of sensors to be a free variable, and enforce a sparse solution, which can be achieved, for example, by $l_1$ regularization technique which yields a convex problem,

\begin{equation}\label{eq:optrelaxl1}
% \nonumber to remove numbering (before each equation)
\min_{\alphav}  J(\G_o(\alphav))+c\|\alphav\|_1  \quad 0\leq \alphav_{tj}\leq 1,
\end{equation}

where the constant $c\geq 0$ is the weighting on the $l_1$-norm penalty. By varying the weight $c$, the number of sensors in the solution set will change to balance the sparsity penalty with the observability measure, tracing the Pareto tradeoff curve between sparsity and observability can be used to achieve a particular level of observability.

§ STRUCTURE GUIDED SENSOR SELECTION

For many systems where the exact network edge structure is unknown or partially known, we can supplement our understanding of the system from structure among the nodes that extends beyond the adjacency structure of the network.
For instance, in constructing social, telecommunications, or postal networks, where the contacts of only some individuals are known, the locations of individuals is informative to network structure <cit.>;
and given our knowledge of the transcription cluster and chromosome territories, the positioning of genes on different chromosomes is important for understanding gene regulatory networks <cit.>.
In each case, we can assign location/positional attributes to each node that form an underlying structure between the nodes, independent of the edge or adjacency structure on the network.

§.§ 1212Small World Networks

networks, characterized by relatively small diameters (distance across the network) and high clustering coefficients (similarity between nodes), are often formed from an underlying structure among the nodes. In 1967, Stanley Milgram's famous experiment, from which the Small World and related six degrees of separation theories stem, mapped the first network by considering the social network between Nebraska and Boston <cit.>; the location of individuals around the United States informed the social network through which packages were mailed to friends and acquaintances. Subsequently, in 1998, Duncan Watts and Steve Strogatz proposed a tunable model, dubbed the Watts-Strogatz model (WS), for generating networks with Small Worlds properties <cit.>. This model introduces long range connections, or “weak ties" to highly structured lattices <cit.>. The initial lattice structure, resembling spatial distribution, serves as the framework for forming Small World networks.

properties have since garnered recognition across a variety of fields where both structural and functional or relational data are captured <cit.>.

Considering the significance of structure in networks and their ubiquity, we study the observability and effective sensor selection methods for such networks. We present results from two simulation to investigate the observability and quality of sensor placement on highly structured, , and random networks, focusing on scenarios where only the lattice or spatial structure is known.

Lattice structures used to create Small World networks are shown with evenly spaced sensor nodes (red) placed across the network. Four types of lattices are shown: (1) ring with 12 vertices, degree of 6, and 3 sensors, (2) $9\times 9$ grid with connections to 4 nearest neighbors and 9 evenly spaced sensors, (3) $9\times 9$ grid with connections to 8 nearest neighbors and 9 evenly spaced sensors, and (4) $9\times 9$ isometric grid with 9 evenly spaced sensors.

§.§.§ 1212Sensor Ranking Distribution from Lattices

Suppose the lattice structure used to generate a network but not the rewired, long-range interactions exist in a network formed from the WS model.

What can the lattice tell us about the expected ranking of each node as a sensor in the network?

To investigate this question, we considered the following model. networks were generated according the the SW model (WattsStrogatz in <ref>) using the lattice structures in <ref>. The unweighted, bidirectional adjacency matrix $\Av$ from each network was used in place of $\Av_t$ in <ref> or <ref>. The value of each node as a sensor was evaluated based upon their contribution to observability as defined by <ref>.

Based upon this model, the expected contribution of each node as a sensor to observability in a small world network mirrors its impact on observability within the lattice structure. In <ref>, the sensor rankings of each node on the lattice are juxtaposed with the aveage sensor ranking over 2500 instances of networks generated from the lattice as well as the sensor rankings on a single network. Qualitatively, the anticipated sensor contributions of each node closely resemble those observed on the lattice.

The rank of each node as a sensor on four different lattice structures and average ranking from SW graphs generated from these lattices is shown. Darker intensity indicates high ranking as a sensors. The top row is the sensor ranking on the lattice shown with black edges; the middle row is the average sensor ranking from 2,500 networks generated from the lattice; the bottom row is the sensor ranking for one instance of a network generated from the lattices.

To asses and provide a scale for the utility of the lattice structure in this sensor selection, we measured the similarity in the sensor selection on the and random networks relative to the lattice structures. As a baseline, we utilized the Erdős-Rényi-Gilbert (ERG) model for random graphs to ensure that the lattice, , and reference point shared identical node and edge counts <cit.>. We recorded the contribution to observability of each node in a vector and calculated the distance between these vectors using the Frobenius norm, as illustrated, as illustrated in <ref>.

The “distance" in sensor rankings between the lattice and networks generated from the lattices are shown. For scale, the distance between the lattice sensor rankings and sensor rankings from random networks are shown. Sensor ranking distance is computed by representing the contribution of each node to observability as entries in a vector, normalizing each network, and computing the distance between vectors from different networks.

Our choice to use the Frobenius norm is arbitrary, so for completion, we evaluate the similarity based upon several additional distances in <ref>-<ref>. Regardless of the choice in norm, the and ERG distributions of similarity in sensor rankings relative to the lattice show clear separation, and across all norms we considered, except for $\|\cdot\|_1$, the sensors were typically more similar to the lattice than random.

Across all norms, the square, diagonal, and isometric lattices typically show more similarity to the sensor ranking on the lattice than the random graph, but this is not true of the ring lattice. In this ring lattice, all nodes are equivalent such that the ring lattice is symmetric up to any rotation of nodes; however, in each of the other three lattices, the structure creates center, edge, and even corner nodes, each with different properties such as their degrees or contributions to observability as shown in <ref>. Similar to the ring lattice, the expected structure of an ERG network is perfectly symmetric. The increased symmetery of the ring lattice explains why across the shown norms, the distributions generated from the ring and random lattices shown in <ref> are similar.

§.§.§ 1212Observability from Spatially Distributed Sensors

Given that the expected sensor ranking of a network can be determined from the lattice, we consider effective stratigies for placing sensors. In line with the supposition that only the lattic structure is known, we consider the peformance of evenly placed sensors on the lattice, as shown in <ref>.

When is placing sensors evenly in space an effective strategy?

To investigate this question, we consider the following experiment. Sensors are uniformly distributed in space on a lattice.

Similar to fig. 1 in <cit.>; however, here we fix $p$ and change the number of iterations whereas Watts' figure was focused on the variability of $p.$

Subsequently, while the sensors remain fixed, we conduct multiple iterations of IterativeWattsStrogatz (see <ref>) wherein edges are progressively rewired, similar to the WS model. Over multiple iterations, the network structure gradually changes from a lattice, to a small world network (see <ref>), and finally a random network, and we evaluate the observability from the fixed sensors nodes throughout. As a baseline for comparison, also we measured the observability of the system given by random sensors.

In the early, highly structured phase, the evenly spaced sensors have increased network observability relative to the randomly spaced sensors, on all latices except for the ring (see <ref>). As the iterations progress and the network approaches a random graph, the utility of the evenly spaced sensors relative to randomly selected sensors diminishes, and the overall observability of the system decreases. The initially increased observability from the structure guided sensors followed by its decrease suggest there exist a set of networks whose sensor selection can be based on the underlying node structure. The underlying lattice structure alone is not sufficient to guarantee good sensor placement, but it does suggest that enforcing constraints to space out the sensors may be beneficial on small world networks.

This figure was made based on the four lattice types with 144 sensors each. Evenly spaced sensors on the lattice improve network observability while the networks are (iterations 0-50), and at later time points, the overall observability and the relative gain in observability from evenly spaced sensors on the lattice as opposed to randomly spaced sensors diminishes .

Generators for Small World Networks

[1]
WattsStrogatz$N, k, p$
$G \gets$ Initialize a ring lattice with $N$ nodes, each connected to its $k$ nearest neighbors
$(u, v) \in$ edges of $G$
random number $< p$
$w \gets$ Randomly select a node (excluding $u$ and its neighbors)
Remove edge $(u, v)$
Add edge $(u, w)$ to $G$
return $G$
IterativeWattsStrogatz$N, k, p, I$
$G \gets$ Initialize a ring lattice with $N$ nodes, each connected to its $k$ nearest neighbors
rewire each edge with probability $p$
return $G$

§.§ 1212Constrained Optimizations

The DSS energy and observability Gramian based sensor selection problems (<ref> and <ref>) can be constrained according to the structure of the system. In the small world example above, we see that placing sensors evenly in space can be an effective strategy for selecting sensors. Here, we modify the DSS optimizations in order to constrain the selection so that at most one sensor per cluster is selected. In principal, these clusters can be created in space or based upon spatial data so that the sensors will be evenly placed in space, but the selection of a clustering method is not of primary importance.

§.§.§ 1212Linear Constraints in Gramian Optimization

Using the trace of the observability Gramian as a measure of observability, the DSS optimization <ref> is formulated as a linear program in eq. (<ref>). In the unconstrained DSS optimization, the only constraint imposed on the linear program is to bound the number of sensor selected. However, we can modify the constraints of this problem based on the spatial clustering of different state variables and solve the optimization problem exactly.

§.§.§ 1212Heuristic Approach for Energy Optimization

To modify the Energy based optimization of <ref>, we use a heuristic, ranking approach. We first solve the maximization problem of <ref> directly and rank the contribution of each sensor at every time. Then, we apply a greedy selection algorithm that selects the top ranked sensor at each time according to their rankings and so as not to violate the rule of selecting multiple sensors per cluster. Pseudocode for this proceedure is provided in <ref>.

Energy based SGSS
[1]
Given: dynamics $\Av_t$ and spatial clusters $C_1,\dots, C_s$
each time point $t$
place genes into a queue according to their contribution to observability from <ref>
the top ranked gene is in a cluster containing no sensors
select the top ranked gene as a sensor at time $t$
remove the top ranked gene from the queue

§.§ 1212Small Worldness of the Genome

Networked models of the nucleus are well established, and the propensity for small world properties have even been seen from molecular dynamics simulations of chromatin conformation <cit.>.
The genome and gene regulatory networks exhibit numerous parallels with classic networks:

* Similar to Milgram's experiment, our current understanding of the full structure in the nucleus and complete set of regulatory interactions of the genome remains incomplete.

* While our understanding of the system is incomplete, like Milgram, recent experimental assays allow us to partially observe both the structure, the location of individuals in Milgram's experiment or the folding of chromatin in the nucleus, and function, social or gene regulatory interactions, of the genome.

* Similar to the WS model, the genome's architecture arises from the polymer structure of chromatin. The modeling and simulation of such polymers employs a bead-on-a-string model that imposes a one dimensional lattice structure while allowing for distant long range interactions or weak ties to form <cit.>.

Despite these similarities, there has been no empirical quantification of the small-worldness of the nucleus from experimental data. Based upon the observability of networks, to apply SGSS for biomarker identification, we examine properties observed in experimental data by (1) constructing networks with adjacency structures that qualitatively resembles Hi-C and (2) quantitatevly evaluating the small world quotient from Hi-C data.

§.§.§ 1212Overview of Experimental Data

Here we provide a concise discussion of the experimental data types we consider in our study.

Hi-C is a genome wide experimental technique to map the spatial arrangement of chromatin within the nucleus <cit.>. Similar to other Chromosome Conformation Capture methodologies, the Hi-C assay (1) cross-links, (2) digests, (3) ligates, and (4) sequences pieces of DNA to identify proximal genomic loci. The resulting data structure is typically viewed as an occurrence or frequency matrix, where the $i,j$th element denotes the number of times the $i$th and $j$th loci of chromatin were observed near one another.

From a data science or controls perspective, Hi-C is both interesting and challenging to consider because it can be analyzed at multiple resolutions, as shown in <ref>. Whereas gene expression occurs in discrete units, Hi-C is analyzed by binning or summing over the base pairs that form the chromatin polymers. Examining genome wide Hi-C at base pair resolution is both challenging, as it requires a matrix with $\approx 9\times 10^{18}$ entries, and unnecessary, as viewing the data at lower resolutions by averaging the number of contacts yields interesting and relevant information. Also, in contrast to eeg or gene expression signals, population Hi-C is less variable throughout time, particularly when the cells are unsynchronized.

Hi-C data can be viewed as a multilevel or multiresolution network. The state variables of Hi-C and the gene regulatory network is not a one-to-one mapping, but there is a mapping between the two.