Information theory for data-driven model reduction in physics and biology

Model reduction is the construction of simple yet predictive descriptions of the dynamics of many-body systems in terms of a few relevant variables. A prerequisite to model reduction is the identification of these relevant variables, a task for which no general method exists. Here, we develop a systematic approach based on the information bottleneck to identify the relevant variables, defined as those most predictive of the future. We elucidate analytically the relation between these relevant variables and the eigenfunctions of the transfer operator describing the dynamics. Further, we show that in the limit of high compression, the relevant variables are directly determined by the slowest-decaying eigenfunctions. Our information-based approach indicates when to optimally stop increasing the complexity of the reduced model. Furthermore, it provides a firm foundation to construct interpretable deep learning tools that perform model reduction. We illustrate how these tools work in practice by considering uncurated videos of atmospheric flows from which our algorithms automatically extract the dominant slow collective variables, as well as experimental videos of cyanobacteria colonies in which we discover an emergent synchronization order parameter. Significance Statement The first step to understand natural phenomena is to intuit which variables best describe them. An ambitious goal of artificial intelligence is to automate this process. Here, we develop a framework to identify these relevant variables directly from complex datasets. Very much like MP3 compression is about retaining information that matters most to the human ear, our approach is about keeping information that matters most to predict the future. We formalize this insight mathematically and systematically answer the question of when to stop increasing the complexity of minimal models. We illustrate how interpretable deep learning tools built on these ideas reveal emergent collective variables in settings ranging from satellite recordings of atmospheric fluid flows to experimental videos of cyanobacteria colonies.

The exhaustive description of a biological or physical system is usually impractical due to the sheer volume of information involved.As an example, the air in your office may be described by a 10 27 -dimensional state vector containing the positions and momenta of every particle in the room.Yet, for most practical purposes, it can be effectively described using only a small number of quantities such as pressure and temperature.Similar reductions can be achieved for systems ranging from diffusing particles to biochemical molecules and complex networks.In all cases, certain relevant variables can be predicted far into the future even though individual degrees of freedom in the system are effectively unpredictable.
The success of these approaches is limited by a fundamental difficulty: before performing any reduction, one has to identify a decomposition of the full system into relevant and irrelevant variables.In the absence of prior knowledge and intuition (e.g. a clear separation of scales), identifying such a decomposition is an open problem [4].It may not even be clear a priori when to stop increasing the complexity of a simplified model or, conversely, when to stop reducing the amount of information needed to represent the dynamical state of a complex system.In both cases one must first determine the minimal number of relevant variables that are needed.The answer to this question depends in fact on how precisely and how far in the future you wish to forecast.Nonetheless, this answer should be compatible with fundamental constraints on forecasting set by external perturbations and finite measurement accuracy [25,26].
In order to address the difficulty identified in the previous paragraph, we develop an information-theoretic framework for model reduction.Very much like MP3 compression is about retaining information that matters most to the human ear [27], model reduction is about keeping information that matters most to predict the future [28,29].Inspired by this simple insight, we formalize model reduction as a lossy compression problem known as the information bottleneck (IB) [30? , 31].This formal step allows us to give a precise answer to the question of how to identify relevant and irrelevant variables.We show how and under what conditions the standard operatortheoretic formalism of dynamical systems [19,32], which underlies most methods of model reduction, naturally emerges from optimal compression.Crucially, our framework systematically answers the question of when to stop increasing the complexity of a minimal model.Further, it provides a firm foundation to address a practical problem: the construction of deep learning tools to perform model reduction that are guaranteed to be interpretable.We illustrate our approach on benchmark dynamical systems and demonstrate that it works even on uncurated datasets, such as satellite movies of atmospheric flows downloaded directly from YouTube and biological datasets composed of microscopy videos of cyanobacteria colonies in which we discover an emergent synchronization order parameter.a State FIG. 1. Interpretable dynamical variables in reduced models via the information bottleneck.(a) The information bottleneck compresses high-dimensional state variables xt, into simpler encoding variables ht providing a controllable trade-off between the degree of compression and the predictive power about the system's future.With deep neural networks, the encoding can be computed directly from data of observed fluid flows (left) or biological datasets, such as fluorescently labeled bacteria colonies (right).In general, the state of the variable xt may comprise time-lagged variables of the intensity field, xt = {It, It+∆t} (right).The amount of compression is determined by the "width" of the bottleneck β [see (1))].The resulting compressed, or encoded, variables ht represent collective variables most predictive of the system's future.(b) Schematic evolution of the encoder p(ht|xt) for varying compression strength β.For low β (high compression), the encoder is trivial and forgets everything about the input xt.After the first IB transition at β1, the encoder becomes non-trivial by gaining some dependence on xt; some features of the input are able to pass through the bottleneck.At subsequent IB transitions, additional features are learned.(c) The point spectrum of the transfer operator contains several slowly decaying modes (red).We show that the most predictive variables that IB systematically extracts correspond to the slowest eigenfunctions of the transfer operator, associated to eigenvalues Λi with |Λi| ≈ 1.In fluid flows, the slowest-decaying eigenfunctions typically represent large-scale coherent patterns of the flow field, while faster-decaying eigenfunctions correspond to variations over shorter length scales.

I. MODEL REDUCTION AS A COMPRESSION PROBLEM
We present here a method to extract collective variables most predictive of the system's future evolution directly from data.This data is composed of a time sequence of measured states x 1 , x 2 , ..., x T .The system state x t could correspond to anything from the position of a single particle to an image of a fluid flow or the fluorescent molecules in a living system (Fig. 1a).The full state can be very high dimensional, with a number of dimensions equal to the number of observed pixels in the case of imaging data.However, the variation of any individual pixel is often of limited interest to us, as noise (either inherent or due to measurement) induces uncertainty about its true value.Individual pixels are, due to this uncertainty, poor predictors of the future state of the system.We can say that they are irrelevant for predicting the future.On the other hand, certain spatially-averaged collective variables may evolve slowly in time, and the future state of the system may be reliably estimated from them.
We seek a way to "encode" each state x t into a simple, lower-dimensional representation h t in a way that isolates these relevant features of the input x t .For instance, in Fig. 1a, both the velocity field of a fluid flow and the images of a dynamic cyanobacteria colony (upper row) may be encoded as a point in a 2D space (lower row).The encoding is given by a probabilistic mapping p(h t |x t ) which can be thought of as a machine which takes a state x t and assigns it to a value h t .The fact that this mapping is probabilistic simply means we may have some uncertainty about the true value of h t even given a measurement of the state.Whether or not some encoding is extracting "relevant" features of the state x t is determined by the extent to which we can use it to predict the future.This predictive power can be quantified by the mutual information between the encoding and the future state, I(H t , X t+∆t ), which tells us how much the knowledge of H t reduces our uncertainty about the future X t+∆t (see Methods).(In our notation, upper-case X t refers to the random variable, while x t refers to a particular value taken by the random variable.)To find a good encoder that extracts relevant features, we might try to find an encoding which maximizes this information.However, an encoder obtained in this way would simply copy the original state, h t = x t , since x t represents all the information we have about the system.In order to encourage the encoder to discard irrelevant features, we simultaneously seek an encoder which maximizes compression by minimizing the information about the original state, I(H t , X t ).This prescription for encoding relevant collective variables can be formalized by the information bottleneck (IB) method [28,30,31].The information bottleneck objective function combines both of our stated goals -compression and predictive power -into one mathematical expression: L IB [p(h t |x t )] = I(X t , H t ) − βI(X t+∆t , H t ). (1) Crucially, the parameter β allows one to tune how much weight to assign to compression versus prediction.For small β the compression term dominates and the optimal encoder is trivial, losing all information about the system.For intermediate β the compression term does not allow X t to be completely captured by H t , so that features of X t must "compete" to pass through to the encoding variable (Fig. 1b).These features are reflected in the form of the encoder p * (h t |x t ) = arg min L IB (2) which provides the optimal trade-off between compression and predictability [29].Our goal is to connect the dynamical properties of the system to the features learned by the encoder.
In any realistic experimental setting, the presence of noise or uncertainty means we cannot predict precisely the future state of a system but instead can only predict a likely distribution of possible future states.Our prediction of the state at ∆t in the future is then represented mathematically as p(x t+∆t |x t ), the probability of observing state x t+∆t given the current state x t .This conditional probability distribution completely characterizes the dynamics of the system, and determines how probability distributions evolve in time: p(x t+∆t ) = p(x t+∆t |x t )p(x t )dx t . ( For Markovian, or "memoryless" dynamics, such an evolution can be understood as the action of a (linear) transfer operator U which acts on probability distributions.U can be decomposed into its right and left eigenvectors as where |ρ n ⟩ are right eigenvectors with eigenvalue Λ n ≡ e λn∆t and ⟨ϕ n | are the corresponding left eigenvectors.λ n are the eigenvalues of the infinitesimal generator of U , known as the Fokker-Planck operator (Fig. 1c).The operator U ess corresponds to the so-called essential spectrum, and we assume that it can be neglected.This is usually possible when the system is subjected to even a small amount of noise, or when some amount of uncertainty is present in the measurements [33,34].The eigenfunctions ϕ n in (4) are in some sense "natural" features of the dynamics, as they evolve independently in time.
Our key observation is that the optimal encoder in (2) can be expressed in terms of the eigenvalues λ n and left eigenfunctions ϕ n of (the generator of) U , where f n (h t ) are factors that do not depend on x t .For an outline of the mathematical steps leading to this see Methods, as well as the SI.These factors effectively determine what the encoder learns about the state x t .In general, there may be a large number of non-zero factors f n so that the learned features are difficult to extract.However, things become simple in the limit of small β, or high compression.When β is small the encoder is trivial: p(h t |x t ) = p(h t ).In this case the value h t is assigned at random with no regard to the state x t of the system.No feature has been learned, and all factors f n are equal to zero.As β is increased, the encoder undergoes a series of transitions at β = β 1 < β 2 < β 3 ... where new features are allowed to pass through the bottleneck (Fig. 1b) [35][36][37][38].
The first transition happens at a finite value of β 1 when the first most predictive feature is learnt.
Surprisingly, we find that at the first IB transition the vector of f n coefficients is dominated by a single term f 1 .
This is our main mathematical result, which we derive by considering a perturbative expansion of the IB objective for small f n .A proof of Eq. 6 with clearly specified technical assumptions may be found in the Methods and SI.
The above statement shows that in the limit of high compression the encoder's dependence on x t is given by the first left eigenfunction ϕ 1 (x t ), which is the slowestvarying function of the state under dynamics given by U .Therefore, Eq. 6 makes precise the intuitive statement that slow features are the most relevant for predicting the future.Our analytical result, while applying only to the dominant eigenfunction, is valid for arbitrary non-Gaussian variables.The question of maximally informative features has additionally been explored in the context of animal vision, where one seeks to understand what features of the field of vision are encoded by retinal neurons [39,40].
We further observe numerically that this picture holds true more generally: also at successive IB transitions, the learned features correspond to successive modes of the transfer operator.This picture is consistent with the exact results known for Gaussian IB, where the encoder learns eigenvectors of a matrix (related to the covariance of the joint X t , X t+∆t distribution) in a step-wise fashion at each IB transition [37].Together, this shows that the most informative features extracted by IB, an agnostic information-theoretic approach, correspond to physically-interpretable quantities -namely transfer operator eigenfunctions.As we show later, the insight above can be leveraged to systematically learn these slow variables directly from data with neural networks [41].

II. INFORMATION DECAY AND THE SPECTRUM OF THE TRANSFER OPERATOR
To develop intuition for information in a dynamical system, we turn to the simple example of a Brownian particle trapped in a confining double-well potential.This might represent, for example, a molecule with a single degree of freedom that transitions between two metastable configurations [42].In the overdamped limit the state of the particle is completely determined by its position X t ∈ R, with dynamics given by the Langevin equation Here, η t is unit-variance white noise, σ controls its strength, and µ controls the shape of the potential V (x).
The deterministic dynamical system undergoes a bifurcation at µ = 0 (Fig. 2a).Sample trajectories, with noise, for a uniform initial distribution of particles are shown in Fig. 2b.For negative µ, the trajectories all converge to a fixed point at x = 0, while for µ > 0 they fluctuate around the fixed points at x = ± √ µ.
To quantify the amount of information about the future state X t+∆t contained in the initial state X t we compute their mutual information (Fig. 2c; see SI for details).The dynamics of X t are Markovian, so that for any sequence of times t 0 < t 1 < t 2 , p(X t2 |X t1 , X t0 ) = p(X t2 |X t1 ).From the data processing inequality, one has [43] I(X t2 , X t0 ) ≤ I(X t1 , X t0 ), which implies that information can only decrease in time.
What governs the rate at which information decays?Here we can already see the role of the spectrum of the dynamics' transfer operator.By exploiting the spectral expansion of the conditional distribution p(x t+∆t |x t ) one finds that for long times the information decays as where expectations are taken over the steady state distribution (see SI). Asymptotically, the information decay is set by the value of λ 1 , the rate of decay of the slowestvarying function ϕ 1 (x) under the dynamics of U .In the limit of infinite time, for any value of µ even weak noise will cause the mutual information to become zero as there is a non-zero probability of hopping between the wells, though this rate of hopping is exponentially small [33].The loss of information in time depends on the bifurcation parameter µ as summarized in Fig. 2d.Note the peak in I(X t , X t+∆t ) for small, positive µ.This corresponds to dynamics where observation of X t strongly informs the future state; recall that the mutual information is maximized when the conditional entropy S(X t+∆t |X t ) ≈ 0 (see Methods).In contrast, for large positive or negative µ, X t is not as informative of X t+∆t even for small times: the initial state is quickly forgotten as the particle approaches the bottom of the single (for µ < 0) or double (for µ > 0) well.
This phenomenon is reminiscent of critical slowing down, which occurs in the noise-free system as µ passes through the bifurcation at µ = 0.For the deterministic dynamics, the slowing down is reflected in the spectrum as a "pile up" of eigenvalues to form a continuous spectrum [33].In the presence of noise, although the continuous spectrum becomes discrete [33,34] there is still a pile-up of eigenvalues characterized by several eigenvalues becoming close to 1 (Fig. 2e).This pile-up gives rise to the information peak seen in Fig. 2d.The peak is not solely due to the closing spectral gap λ 1 − λ 2 , but is also impacted by the subdominant eigenvalues which accumulate at µ ≈ 0.2 (SI Fig. S4).

III. KNOWING WHEN TO STOP
For discrete encoding variables h, the information bottleneck partitions state space and reduces the dynamics on x to a discrete dynamics on h.Such reductions of complex systems to symbolic sequences via partitioning of state space has attracted attention for more than half a century in both theoretical and data-driven contexts [44][45][46][47][48][49][50].Several works have approached this partition problem from a dynamical systems perspective, linking "optimal" partitions to eigenfunctions of the (adjoint) transfer operator [25,51].In this setting, a central question is "when to stop" [25,26,48,49]: how many states does h need in order to capture statistical properties of the original dynamics?
We consider this question by finding the optimal IB encoder in the limit of low compression, β ≫ 1, but fixed encoding capacity N H (where H t ∈ {0, ..., N H − 1}), i.e. the encoder is only restricted by the number of symbols it can use.An analogous setup was used in the context of renormalization group (RG) transformations in [52][53][54], which results in effective model reduction due to the "sloppiness", or irrelevance, of certain system variables [55,56].In this regime, the encoder learned by IB is deterministic; we are learning an optimal hard partition of state space.This can be seen by noting that I(H t ; X t+∆t ) = S(H t ) − S(H t |X t+∆t ) is maximized when the latter term is zero, which happens when x t unambiguously determines h t , i.e. when p(h t |x t ) ∈ {0, 1} for all x t .The details of how the encoder is computed are discussed in the next section.Fig. 2f shows that the number of states necessary to describe the system depends strongly on the value of µ.For |µ| ≫ 0, a two-state discrete variable h t ∈ {0, 1} suffices to describe the system's future.Increasing the number of reduced variables N H does not allow more information to be captured.Near the information peak at µ ≈ 0.2 this changes: predicting the future state of the system requires a more complex hidden variable of up to N H ≈ 10 values.Above, we saw that this peak arises due to the pile-up of eigenvalues at µ ≈ 0.2.The content of the transfer operator spectrum is thus directly reflected in the number of encoding variables needed to capture the system's statistics.
Noise can have a similarly dramatic impact on the reduced model complexity.Indeed, noise in some form, either inherent to the dynamics or due to measurement error, is necessary for a model to be reducible.In purely deterministic systems where the future state is a bijective function of the present state, information does not decay and complete knowledge of the state is required to predict the future.
Consider a fluctuating Brownian particle as above, where now each of the wells is split into two smaller wells, giving a total of four potential minima (Fig. 3a).As the system is in steady state, the standard deviation of the fluctuations σ corresponds to an energy scale E σ = σ 2 = 2k B T .For small E σ , the system rarely transitions between the four potential minima.In this case, knowledge of the initial minimum is very informative of the future state of the particle.In contrast, for large fluctuations the particle can spontaneously jump between shallow minima in each large well, so that the system immediately forgets about the precise potential minimum it was in.Information about the shallow minima has been "washed out", and only the information about the larger double-well structure remains.
To see this reflected in the information, we again consider an encoding of the initial state into a discrete variable H t ∈ {0, ..., N H − 1}.In both the small and large noise scenarios, a variable with N H = 2 encodes approximately one bit of information (Fig. 3b), corresponding to an H t which distinguishes the two large wells for x ≶ 0. For large noise this is essentially all the information that can be learned; increasing the capacity of the encoding variable beyond this provides only marginally more information about the future state (Fig. 3c).In the small noise case, the information between the encoding and the future state continues to increase to approximately two bits at N H = 4, after which it plateaus.The encoding has learned to distinguish each of the four potential wells.These observations are reflected in the transfer operator spectrum shown in Fig. 3d.For small noise, the eigenvalue λ = 1 is four-fold degenerate which indicates the existence of four regions that can evolve independently under U , giving rise to four steady state distributions satisfying U ρ = ρ.These regions correspond to the potential minima.Hops between the separate minima are exceedingly rare, so that the dynamics essentially take place in the four minima independently.With larger σ the degeneracy is lifted, resulting in one dominant subleading eigenvalue followed by a gap.The corresponding eigenfunction is one which is positive (negative) on the right (left) side of the large potential barrier at x = 0: the only relevant piece of information is which of the large wells the initial condition is contained in, and all other information is lost exponentially quickly.

IV. TRANSFER OPERATOR EIGENFUNCTIONS ARE MOST INFORMATIVE FEATURES
Until now we have concerned ourselves with encodings whose capacity is limited only by the number of variables, rather than by the compression imposed by a small value of β.In the regime of small β, or high compression, features of the state x t are forced to compete to make it through the bottleneck h t .By studying the behavior of the encoder in this regime, in particular its dependence on x t , we may identify the most relevant features of the state variable and show that they coincide with left eigenfunctions of the transfer operator.
We return to the simple example of a particle in a double well with dynamics given by (7) which we map to a discrete variable H t ∈ {0, ..., N H − 1}.In this system the IB loss function (1) can be optimized directly, as shown in Ref. [30], using an iterative scheme known as the Blahut-Arimoto algorithm [43] (see SI).
To focus on the properties of encodings for varying degrees of compression β, we consider a fixed set of dynamical parameters µ and σ.Increasing β reduces the amount of compression, i.e. "widens" the bottleneck, allowing more information to pass into the encoder.This leads to a series of IB transitions which are sketched in Fig. 1b and shown quantitatively in Fig. 4a.The form of the optimal encoder changes qualitatively at these transitions.Before β 1 , the optimal encoder has no dependence on x so that p(H t = h i |x t ) = const for all h i .After the first transition, the encoder begins to associate regions of x to particular values of h.We are interested in the form of the encoder at β ≳ β 1 , just above the first IB transition, as this reflects the most informative features of the full state variable x (Fig. 4a).The dependence of p(h t |x t ) on x can be explained by a stability analysis of the IB Lagrangian (see Methods and SI).Stability is governed by the eigenvalues η i of the Hessian of the IB Lagrangian with respect to the parameters f n (h t ) in Eq. 5.These parameters tell us how much the encoder "weights" each transfer operator eigenfunction; f n (h t ) = 0 (for n > 0) corresponds to the uniform, or trivial encoder p(h t |x t ) = p(h t ).
For small β all eigenvalues η i are positive, indicating that the uniform encoder is a stable minimum of the IB Lagrangian.In Fig. 4b we show the smallest two eigenvalues of the IB Hessian when evaluated at the uniform encoder.At the first transition one eigenvalue becomes negative, so that the uniform encoder is unstable.The eigenvector corresponding to the unstable eigenvalue η 1 indicates how the weights f n (h t ) should be adjusted to lower the value of the IB Lagrangian.Our numerics confirm that these weights are dominated by f 1 as expected from our analytical result (Fig. 4c, top).By taking the logarithm of the encoder after the transition, we can in- dependently confirm that the encoder depends only on ϕ 1 (x) (Fig. 4d).
Our stability analysis predicts that a second mode becomes unstable at the second IB transition β ≈ β 2 (Fig. 4b).Here we see that this unstable mode selects f 2 , and that the encoder correspondingly gains dependence on ϕ 2 (x) (Fig. 4d).Note that in general, η 2 must not necessarily become negative precisely at β 2 because the stability analysis is performed at the uniform encoder while the true optimal encoder has already deviated from uniformity.In the SI, we perform the same analysis for a triple-well potential where this difference is more apparent.

V. DATA-DRIVEN DISCOVERY OF SLOW VARIABLES
IB finds transfer operator eigenfunctions by optimizing an information theoretic-objective that makes no reference to physics or dynamics.This suggests it may be used for the discovery of slow variables in situations where one lacks physical intuition.The utility of exact IB for this purpose is limited because it requires knowledge of the exact conditional distribution p(x t+∆t |x t ) which is difficult to estimate in many real-world scenarios.Fortunately however, the IB optimization problem can be replaced by an approximate variational objective introduced in Ref. [41] that can be solved with neural networks.We refer to this as variational IB.In the remainder of this paper, we show how to implement these networks for the discovery of slow variables directly from data.
First we show numerically that the results of the previous sections remain valid for high-dimensional systems by considering a simulated data of fluid flow past a disk [57].The state of the system is given by a two-dimensional velocity field v(x) ∈ R 2×N pixels , where N pixels ∼ O(10 5 ) (Fig. 5a).Fluid flows in from the left boundary with a constant velocity v 0 êx past a disk of unit diameter.At Reynolds number Re ≳ 150, the fluid undergoes periodic vortex shedding behind the disk, forming what is known as a von Kármán street.
What do the true eigenfunctions look like in this system?Because it is well approximated by linear dynamics, eigenfunctions of the adjoint transfer operator are linear functions of the state variable, where m (n) is the n-th mode (often referred to as a Koopman mode [19]) and angled brackets denote integration over space.The true eigenfunction and corresponding modes can be computed via dynamic mode decomposition (DMD) [11,12], as described in the SI.The eigenfunctions for this system are in general complex, and come in conjugate pairs: ϕ 2 (x) = ϕ * 1 (x).In this situation any linear combination of ϕ 1 and ϕ 2 will decay at the same rate, and hence we expect to learn some arbitrary combination of the two dominant eigenfunctions, or equivalently a combination of the real and imaginary parts of ϕ 1 .We therefore take a two-dimensional encoding variable [h 0 , h 1 ], so that it can represent the full complex eigenfunction rather than only the real or imaginary part.
Our learned latent variables are oscillatory with the correct frequency as shown in Fig. 5b-c.A more stringent test is whether we are also learning the correct mode m (1) .From the true eigenfunctions, the modes can be extracted by computing the gradient As the learned function h[v] is a neural network, we can efficiently compute gradients of the network with respect ) and from VIB (m (IB) ).Koopman modes from VIB are computed as gradients of the latent encoding variables as described in the main text.Red corresponds to positive values and blue to negative; the magnitudes of the modes are not directly comparable.
to the input field where we have separated the part of the gradient which is independent of v from a residual part which is dependent on v.If h corresponds to the true eigenfunction, we expect that m (IB) is approximately equal to the Koopman mode m (1) , and that g res is small.We indeed find this to be the case; Fig. 5d shows these gradients averaged over several instantiations of the neural network, which corresponds strongly to the true mode.Details concerning both the averaging procedure and the residuals g res can be found in the SI.This shows that variational IB not only recovers the essential oscillatory nature of the dynamics, but does so by learning the correct slowly varying functions of the state variable given by the adjoint transfer operator eigenfunctions.

VI. RELEVANT VARIABLE IDENTIFICATION IN LABORATORY-GENERATED AND ATMOSPHERIC FLOWS
The scenario above is characterized by high-dimensional data and few samples; training was performed with only ∼ 400 samples.We now demonstrate that our framework continues to hold approximately and yield interpretable latent spaces even for real-world fluid flow datasets scraped directly from videos on Youtube [58,59] (Supplementary Movie 1).
The first shows a von Kármán street which forms as water passes by a cylindrical obstacle at Reynolds number 171, with flow visualized by a dye injected at the site of the obstacle [58].We take a background-subtracted grayscale image of the flow field as our input (Fig. 5e) and task VIB with learning a two-dimensional latent variable as above.Also here, variational IB learns oscillatory dynamics of the latent variables (Fig. 5f).We visualize the function learned by the encoder by considering gradients of the latent variables, which show the same structure as those obtained for the x component of the simulated data (Fig. 5g).This is expected, as the x-component of the velocity field has similar glide reflection symmetry as the intensity image.
Next, we apply variational IB to a von Kármán street arising due to flow around Guadalupe Island, which was imaged by a National Oceanic and Atmospheric Administration (NOAA) satellite [59] (Fig. 5h).The video consists of only 62 frames, and the von Kármán street undergoes a single oscillation.Even with this small amount of data, the variational IB neural network learns latent variables which capture this oscillation and have the expected dependence on the input variables (Fig. 5i-j).As in the first experimental example, the gradients of the encoding variables show the glide symmetry of m x due to the symmetry of the intensity pattern in Fig. 5d.This symmetry is less clear in the component ∂h1 ∂I(x) , which is likely due to the fact that the von Kármán street is not as fully formed in this data as in our previous examples.

VII. RELEVANT VARIABLE DISCOVERY IN CYANOBACTERIAL POPULATIONS
We now demonstrate how variational IB may be used as an aid for collective variable discovery in situations where physical intuition may not be a useful guide -collective behavior of biological organisms (Supplementary Movie 2).Here, we ask what the most predictive variables are for predicting the evolution of populations of cyanobacteria (Synechococcus elongatus).The dynamics of the colonies are driven by several factors: growth and division of individual bacteria, translational motion of groups of bacteria as they are pushed by their neighbors, as well as the circadian oscillations within each bacterium (Fig. 6a).These oscillations are controlled by three Kai proteins [60] and depend in particular on the ratios of the copy number of these proteins which can be tuned experimentally [61].
We were provided with videos of 10 cyanobacteria colonies that were grown under various conditions that impact their dynamics.However, as a test of our method, we were blinded to these conditions until we had performed our analysis.The videos are sequences of fluorescent images, taken once per hour, which show the clock state of each individual bacteria visualized with a fluorescent marker EYFP driven by the kaiBC promoter.Here, we focus on collective variables which are predictive of the state of the interior of the colony and not the growth in area of the colonies.We therefore crop the images to the interiors of each colony (SI Fig. S9).This allows us to isolate the motion of individual bacteria and fluorescence oscillations (Fig. 6b).
Our input to the variational IB neural network are these cropped images augmented with a time-lagged image of the same region (Fig. 6c).The purpose of this time lag is to make the dynamics Markovian: due to the oscillatory intensity field, if one observes only a single time point it is unclear whether the intensity is currently increasing or decreasing.These time lagged pairs comprise our system state, X t = {I(x, t), I(x, t + τ )}, where τ is the duration of the time lag.Here we take τ = 3 hr and a prediction time horizon ∆t = 8 hr, but find that choosing different ∆t or τ does not change our results (SI Fig. S9).
With variational IB we compress the state X t into a latent variable h of variable dimension (Fig. 6d-f).We train the neural network on the entire dataset of all 10 colonies simultaneously.The dynamics in latent space undergo clear oscillations, indicating that the relevant variables encode primarily the intensity fluctuations rather than, for example, the spatial locations of the bacteria.Notably, the trajectories are essentially two-dimensional, even when the encoding space is higher dimensional.This is reflected in the information retained about the future state, I(X t+∆t , H t ).We see that increasing the dimension of the embedding space beyond two leads only to marginal increases in I(X t+∆t , H t ); this tells us "when to stop" (Fig. 6f).We independently verify this by using principal component analysis to characterize the geometry of embedded trajectories, and find that even in higher dimensions the trajectories occupy a two dimensional subspace (SI Fig. S10).In the following, we therefore restrict our focus to the dim H = 2 case.
We noticed that there were notable differences in the radius of latent space oscillations from colony to colony, two of which are highlighted in (Fig. 6g).To understand this difference, we examined the original microscopy time series corresponding to both large and small latent radius (Fig. 6k) and found that while the large-radius sample showed clear, nearly uniform oscillations in intensity, the small-radius samples appeared much more heterogeneous.
To quantify this we consider each pixel to be an independent oscillator, akin to a spatial Kuramoto model [62][63][64], and compute a global synchronization order parameter r(t) (see SI).For each colony we calculate the timeaveraged synchronization ⟨r(t)⟩ t and find that two clusters emerge corresponding to high and low synchronization (Fig. 6i).These clusters are precisely those representing trajectories of large and small latent radius (Fig. 6j), suggesting that variational IB learns to encode the synchronization of the colony in the latent variable radius.As a check, we perform IB on a simulated locally-coupled Kuramoto model as a system which shares many features of the experimental system.Here we also learn an encoding in which the latent radius corresponds to the synchronization order parameter (SI Fig. S11).In the SI we compare the performance of variational IB to several other model reduction methods and find that IB delivers more interpretable and well-behaved features.This is likely due to the fact that many standard methods for data-driven model reduction rely on assumptions about the dynamics which may not be appropriate in the case at hand, such as linearity.Even among deep learning methods such as time-lagged autoencoders that are free of such assumptions, the variables learned by IB appear more interpretable.This increased interpretability is likely due to the compression term which effectively regularizes the latent space by encouraging the network to learn slow transfer operator eigenfunctions.While there are many specific use variants of DMD [13,[65][66][67][68] or autoencoders for dynamics [20][21][22]69] that may outperform variational IB in some cases, we find that in this real-world example it yields the smoothest and most interpretable latent variables without requiring tailored pre-processing steps (SI Fig. S12).By using variational IB we could reduce a complex system with multiple dynamical components -cell growth, division, and gene expression fluctuations -into a low dimensional form that retains only the most relevant information for the future.In addition to the insight that the dynamics are dominated by oscillations in two dimensions, the latent variables clearly distinguished trajectories into two groups that were not apparent a priori.We were provided this data as a "blind" test with no knowledge of the underlying system.After we performed our analysis, it was revealed to us that these bacterial colonies have been engineered to control the translational efficiency of the Kai proteins by varying theophylline concentration [61].The synchronization order parameter discovered by variational IB corresponds to differing experimental concentrations of theophylline, which is in agreement with the findings in Ref. [61].IB can thus serve as a way to connect experimental control parameters to effective changes in dynamics.

VIII. CONCLUSION
We have related information-theoretical properties of dynamical systems to the spectrum of the transfer operator.We illustrate our findings on several simple and analytically tractable systems, and turn them into a practical tool using variational IB, which learns an encoding variable with a neural network.The latent variables of these networks can be interpreted as transfer operator eigenfunctions even though the network was not explicitly constructed to learn these: it optimizes a purely information-theoretic objective that contains no knowledge of a transfer operator or dynamics.This allows one to harness the power of neural networks to learn physicallyrelevant latent variables.Biological systems are an ideal setting for such methods: despite their apparent complexity, they can often be captured by low-dimensional descriptions which are difficult to identify by physical considerations alone [47,[70][71][72].We have shown that variational IB is a potentially powerful tool for these cases, and can discover slow variables even directly from image data without significant preprocessing.

Mutual information and entropy
Let X be a random variable which takes values x that are observed with probability p(x).The entropy of this distribution measures the predictability of the outcome of a measurement of X and is given by [43] S(X) = − dx p(x) log p(x).
Given another random variable Y , such that X and Y have a joint distribution p(x, y), we can ask how much information is shared between these two variables.This is given by the mutual information This can be interpreted as quantifying how much (on average) a measurement of Y can reduce our uncertainty about the value of X (( 12)).

The information bottleneck
The information bottleneck [30] is an example of a rate-distortion problem which seeks to find an optimal compression which minimizes some distortion measure with the original signal [43].Concretely, we call X the source signal, and let H denote the compressed signal.In IB, rather than using an a priori unknown distortion function, one seeks to ensure that the compression retains information about an additional relevance variable Y .As noted in the main text, the IB optimization objective is given by the Lagrangian where in our case the source signal X is the state of the system X t at time t, and the relevance variable is the state of the system X t+∆t at a future time t + ∆t.The encoder which optimizes this objective can be solved for exactly and is given by [30] p Encoder in terms of transfer operator eigenfunctions To connect the optimal encoder to the transfer operator, we first rewrite (15) in terms of the transition probabilities, where we have absorbed terms in the exponent which depend only on h t or x t into the normalization factors.Into the above equation, we replace the transition probability with the spectral decomposition From this, the (5) of the main text immediately follows, where which may be interpreted as a sort of cross entropy (ρ n is generally not a probability distribution) between each right eigenfunction and the decoding of h t into the future state x t+∆t .
To study the behavior of the encoder in the limit of high compression, we consider a transfer operator U with infinitesimal generator L U .For L U with a discrete spectrum with eigenvalues satisfying 0 = λ 0 > λ 1 > λ 2 ≫ λ 3 ... and for β just above the first IB transition β 1 , we show that the optimal encoder is given approximately by with corrections due to the second eigenfunction given by ) where Γ = λ 1 − λ 2 > 0 denotes the spectral gap.To see this, we compute the Hessian of the IB Lagrangian Here we assume a finite alphabet of size N H , i.e. h µ with µ ∈ {1, ..., N H }. At the uniform encoder, i.e. f n (h) = 0, the Hessian decomposes into a tensor product where A β depends only on the indices of the coefficients and G captures the dependence on h ν .The only part which depends on β is A β .We are concerned with the sign of the eigenvalues of H.A negative eigenvalue indicates that L IB is unstable to a perturbation in f , which means the loss can be lowered by changing f away from the trivial encoder at f n = 0.Because the eigenvalues of a tensor product of matrices are products of the eigenvalues of the component matrices, the eigenvalues of H change sign when those of A β do.A β and its spectrum can be computed, which we do in the SI.The result of this calculation is that the first eigenvalue to become negative is associated with the eigenvector v = (1, 0, 0, ...).This computation is exact for equilibrium systems, which are those in which the steadystate flux vanishes, but in nonequilibrium systems there may generally be a correction proportional to the flux.In summary, this means that in the limit of high-compression only the first component f 1 becomes non-zero, hence the encoder has the form given in (19).

Variational IB compared to other dimensionality reduction techniques
Variational IB (VIB) is by no means the only numerical method for performing data-driven model reduction.Here we provide a brief overview of the benefits and shortcomings of VIB with respect to other methods; an extended discussion can be found in the SI.
One class of methods is based on linear projections, such as principal component analysis (PCA), dynamic mode decomposition (DMD) [11,12], or (time-lagged) independent component analysis (TICA) [10] (which is equivalent to DMD [73]).These methods can be extended to take into account non-linearity by introducing a library of non-linear terms on which one then applies the above methods, such as in kernel PCA [74] or extended DMD (eDMD) [13].These methods have the advantages, relative to VIB, that their optimization (even for the extended algorithms) relies only on linear projections which are fast and interpretable.However, the success of these methods depends on the choice of an appropriate library of functions so that the projection onto this space is closed under the dynamics.Choosing an appropriate library is not always possible [75,76].
A second category of non-linear dimensionality reduction techniques are graph-based or similarity-based methods, which typically assume that the data is distributed on a low-dimensional manifold embedded in a higherdimensional space [77,78].One prominent example is diffusion maps [14], which starts from a set of data snapshots and, assuming the system evolves diffusively on short times, constructs an approximate transition matrix from which one can compute eigenfunctions to parameterize the data manifold.The assumption of diffusive dynamics can be violated when data is not sampled sufficiently frequently.This likely explains our finding that VIB produced more well-behaved low-dimensional embeddings on the cyanobacteria dataset (SI Fig. S12).VIB has the additional advantage, relative to this and similar methods, that it explicitly takes dynamics into account without the strong assumptions required by diffusion maps.
Finally, deep neural networks can be used for model reduction through encoder-decoder architectures that attempt to reconstruct the data from a low-dimensional latent space; VIB falls into this class of methods.Some standard neural network architectures from this class include autoencoders (AEs) and variational autoencoders (VAEs).For dynamical systems in particular, extensions to these methods have been proposed which impose constraints on the latent dynamics, such as linearity [20][21][22]69].Autoencoders often produce poorly-behaved latent spaces that distribute the latent variables on a narrow manifold with sharp features, see for example [69].By regularizing the latent embedding to encourage smoothness, variational autoencoders can remedy some of these issues.We note that the VIB loss is very similar to a VAE loss with the contrastive InfoNCE loss replacing the reconstruction loss, so we expect that for many problems these should perform similarly.Other dynamically-constrained architectures such as in [20][21][22] work well for deterministic systems but it is unclear what effect stochasticity has on their performance.In our examples we have seen that VIB works well on noisy data.
In general when investigating a new system it is good practice to start by attempting to perform dimensionality reduction with linear methods such as PCA or DMD because they are fast, straightforward to implement, and easy to interpret.In situations where linear techniques are not sufficient, VIB may be preferable to other methods because it is guaranteed to find dynamically relevant variables (in contrast to diffusion maps, t-SNE, AEs, VAEs, etc.) and it does not require that one performs the carefully tailored preprocessing steps that are required by eDMD or kernel PCA, or other variants of DMD [65][66][67][68].Additionally, it works well even when the dynamics are highly stochastic as shown in the cyanobacteria dataset.
Appendix A: Information bottleneck The information bottleneck was originally formulated in Ref. [30] as a rate distortion problem.Rate distortion theory describes how to maximally compress a signal (i.e.minimize its communication rate) such that it remains minimally distorted, which, however, presupposes the knowledge of a distortion function specifying which features need to be preserved [43].In Ref. [30], Tishby et al. introduce a variant of the problem where the a priori unknown distortion function is not used, but rather one seeks to ensure that the compression retains information about an auxiliary variable Y correlated with the signal.As the correlations with Y define implicitly the relevant features of the signal to be preserved, Y is called the relevance variable.Concretely, we call X the source signal, and H denotes the compressed signal.The random variables form a Markov chain H ↔ X ↔ Y , meaning that H and Y are conditionally independent given X: p(y|h, x) = p(y|x)p(h, x) p(h|x, y) = p(h|x)p(y, x) As noted in the main text, the IB optimization objective is given by the Lagrangian To enforce normalization of the p(h|x) one introduces a Lagrange multiplier λ(x) so that the full optimization function is The encoder which optimizes this objective can be solved for exactly.Using the following functional derivatives, one can compute the derivative of the Lagrangian.By evaluating the derivative and setting to zero, one finds which can be rearranged to give the optimal encoder By absorbing terms which only depend on h and x into p(h) and N (x), respectively, the encoder can be expressed as When β < 1, it follows from the data processing inequality I(X, H) ≥ I(Y, H) that (A1) is minimized by a trivial encoder p(h|x) = p(h).In this case, L IB = 0 and no information passes through the bottleneck.As β is increased, more information is allowed through the bottleneck until the relevant variables begin to gain a dependence on the state x.This occurs suddenly for a certain value of β = β 1 > 1 at which the encoder becomes non-uniform and I(H, X) becomes non-zero.This is referred to as an IB transition; for increasing β, there may be a sequence of transitions at β 2 , β 3 , ... etc.

Appendix B: Numerically solving "exact" IB and Ulam approximations of the transfer operator
The optimal IB encoder (A3) can be found using the Blahut-Arimoto (BA) algorithm [43].As described in detail in Refs.[30,79] and sketched in Fig. S1, the algorithm is an iterative procedure, where the encoder at iteration k + 1 is Blahut-Arimoto Compression Info.max FIG.S1.Exact IB.Exact IB finds an optimal encoder using an iterative Blahut-Arimoto algorithm, described in the main text, for a rate-distortion problem with Kullback-Leibler distortion [30].This process requires access to (an estimate of) the transfer matrix and the steady state distribution.Once an optimal encoder has been found, relevant quantities such as mutual information can be computed.

updated according to
The first line simply plugs the previous estimate of the encoder in to Eq. D2 (and normalizes the distribution), while the following two lines update marginal and conditional distributions using the new estimate of the encoder.In Refs.[30,79] it is shown that this algorithm converges.
The BA algorithm requires access to the conditional distribution p(x t+∆t |x t ) for each x t , as well as the steady state p(x t ).To solve the IB optimization problem we therefore need a numerical approximation of the transfer operator, which we obtain by an Ulam approximation [44,45].In brief, one divides space into bins and computes a finite-dimensional approximation to the conditional distribution as where N i is the number of trajectories starting in bin i and N i→j the number of observed transitions from bin i to bin j.The transfer operator is then approximately given by and eigenvalues and eigenvectors of U can be computed by diagonalizing P .

Appendix C: Transfer operators
In this work we consider dynamical systems given by a Langevin equation where the first term is the deterministic part of the dynamics, and the second term corresponds to noise, where This framework captures purely deterministic dynamics, which are obtained by setting D = 0.The transfer operator describes the evolution of probability distributions.For dynamics given by a Langevin equation, probability distributions evolve according to the corresponding Fokker-Planck equation Here we recognize L as the infinitesimal generator of the transfer operator U : evolution for a time t is given by U = e tL .
The adjoint operator is given by the so-called backward Kolmogorov equation, and it describes the evolution of functions ϕ.
For D = 0, L is the generator of the Perron-Frobenius operator, while its adjoint is the Koopman operator [19].

Appendix D: Derivations
Here we derive the form of the optimal encoder at the first IB transition, In Ref. [30], it is shown that the encoding minimizing the IB objective satisfies the implicit equation where N is a normalization factor ensuring dh p * β (h|x) = 1, and The Kullback-Leibler divergence between two probability distributions p(x) and q(x) is defined by in the D KL in (D2), the integration is done with respect to the common random variable X t+∆t of the two conditional distributions.Equation D2 contains the conditional distribution p Xt+∆t|Ht , which is the composition of two operations: decoding (going from H t to X t ) and time evolution (going from X t to X t+∆t ).The decoding is performed using Bayes' theorem p Xt|Ht p Ht = p Ht|Xt p Xt , so that we have Note the appearance of the encoder p Ht|Xt , making (D2) only an implicit solution of the optimization problem.The conditional distribution p Xt+∆t|Xt can be written in terms of right and left eigenvectors of U .Here we neglect U ess , which is the operator corresponding to the essential spectrum.The essential spectrum is the part of the spectrum that is not the discrete spectrum; by neglecting it, we are assuming that the essential radius ρ ess (the maximum absolute value of eigenvalues in the essential spectrum) is small enough compared to the first few eigenvalues |λ n | in the point spectrum.While purely deterministic systems may exhibit a large essential radius, the introduction of noise causes the essential spectrum to shrink or disappear [33,34].We can then write the conditional distribution in terms of right and left eigenvectors of U , The Kullback-Leibler divergence in (D2) can be represented as a sum of two terms, The first term has no dependence on h and can hence be absorbed into the normalization N (x).Plugging the decomposition (D6) into the second term leads directly to where which can be understood as a quasi-cross entropy between the right eigenvector ρ n and the distribution on x t+∆t obtained by decoding h.Equation D8 depends only on the existence of the spectral decomposition; we have made no other assumptions on the dynamics (beyond Markovianity) up until this point.

Perturbation theory at the uniform encoder
In this section we derive our main result for the optimal encoder in the limit of high compression (small β).
Mathematical result -Consider a transfer operator U with infinitesimal generator L U .Assume that L U has a discrete spectrum with eigenvalues satisfying 0 = Re λ 0 > Re λ 1 > Re λ 2 ≫ Re λ 3 ....Then, for β just above the first IB transition β 1 so that f n (h) → 0, the optimal encoder is given approximately by with corrections due to the second eigenfunction given by f 2 (h) ≈ f 1 (h)e −Γ∆t + O(e −2Γ∆t ) where Γ = λ 1 − λ 2 > 0 denotes the spectral gap.
To show this we compute the f n (h) coefficients at the onset of instability.The instability at the first IB transition at β = β 1 corresponds to the emergence of negative eigenvalues in the Hessian of the IB loss, Similar perturbative approaches have been studied in other contexts within IB in [35,38,80].For concreteness, we consider a discrete state x ∈ {1, ..., N X } and discrete encoding h ∈ {1, ..., N H }. We express the encoder as in Eq.D8 where we now expand the normalization terms In the following, we simplify notation as where φn (x i ) = βe λn∆t ϕ n (x i ).For the uniform encoder, f µ n = 0 for all n ̸ = 0. Hereafter we assign φ0 = 1 f µ 0 , taking care to track the correct Greek index which otherwise is not present for φn (in the following, n is always paired with µ and m always with ν).For example, a single derivative of the marginal is given by To compute the Hessian, we use the fact that λ ∂ µ n p(h λ ) = 0 (by exchanging sum and derivative), and that p(h µ |x i )| f =0 = p(h µ ).We note in passing that we find the first derivative terms of L IB vanish, see also [36,38].The compression term is given by and the information maximization term is Angled brackets with no subscript correspond to an average with respect to the steady state distribution, ⟨•⟩ = xi • p(x i ).For the second term, we retain the integration variable y for clarity.In sum, the Hessian of the Lagrangian is given by which lives in R (N X ×N H )×(N X ×N H ) ; we index H with a multi-index (µ, n).From the form of Eq.D15, we see that H is given by a Kronecker (tensor) product Eigenvalues of H are given by products of eigenvalues of A β and G, while eigenvectors are given by the tensor product of eigenvectors of A β and G.The appearance of unstable directions of the Hessian corresponds to the appearance of negative eigenvalues in its spectrum.As G does not depend on β, zero-crossings of eigenvalues of H therefore correspond to zero-crossings of eigenvalues in the spectrum of A β .
Equilibrium systems: We first study the stability of the matrix A β nm in the case of quasi-equilibrium dynamics.In particular, we consider a Fokker-Planck operator L FP of the form The steady state distribution ρ 0 (x) satisfies where J i = f i ρ 0 − D∂ i ρ 0 is a flux.Our main assumption in this subsection is that in the steady state, fluxes vanish: J i = 0.In this case, left eigenfunctions ϕ n of L FP become right eigenfunctions when multiplied by ρ 0 .To see this, note that where the first term disappears because ρ 0 is the steady-state distribution, and the third term disappears because of our above assumption on disappearing fluxes, J i = 0.The following inner product then satisfies This shows that if λ n ̸ = λ m , the inner product vanishes.Similarly, one can show that the same must be true for ⟨ ρn ρ0 , ρm ρ0 ⟩.These are precisely the types of terms appearing in the Hessian.Consequently, for equilibrium (no flux) systems the Hessian is diagonal.The first term follows directly from the above, while the second is given by The full matrix A appearing in the Hessian then takes the form except for the n = m = 0 term, which is A 00 = 0.The eigenvalues are given directly by these diagonal elements.For small β these are all positive, and become unstable one after the other at which are increasing in order (remember λ i ≤ 0).It follows that at the first transition, only f 1 becomes non-zero, and hence the encoder takes the form given by (D1).In the equilibrium case, the encoder learns exclusively the first eigenfunction, with no correction due to the second eigenfunction.General Case: We now show that the first component f 1 is selected at the first IB transition even when the flux J is non-zero.From our assumption on the spectrum of L U , namely 0 > Reλ 1 > Reλ 2 ≫ Reλ 3 ... it follows that near the first IB transition the IB loss is given by with a 12 = e (λ1+λ2)∆t ⟨ϕ 1 ϕ 2 ⟩ − βB 12 = e (λ1+λ2)∆t â12 where we have introduced the shorthand and all angled brackets denote averaging with respect to the steady state distribution.The stability of the uniform encoder at f n = 0 is given by the stability of the 2×2 matrix above.The eigenvalues ω i and eigenvectors v i of this matrix can be computed explicitly, In what follows, we will express quantities in terms of the spectral gap Γ = λ 1 − λ 2 > 0. For example, the expression above can be written The first eigenvalue to become negative is η + .We are interested in the ratio of the corresponding eigenvector components, , which tells us how much the encoder will depend on the first eigenfunction ϕ 1 (x) compared to the second.This ratio is given by This suggests that we must make one additional assumption, which is that the factor â11 /â 12 is not small.This is true whenever the flux is small, as Similar to the equilibrium case, we see that at the first transition the encoder depends only on f 1 , giving (D1).In contrast to the equilibrium case, there may be a small correction due to the second eigenfunction, however this becomes exponentially small for long times ∆t.
Note that in this calculation, ϕ 0 (x) is constant (following from the assumption of a non-degenerate eigenvalue at 0) so that the f 0 (h) factors can be absorbed into p(h).This changes in the case of a degenerate ground state, which corresponds a situation where there are decoupled sectors in which the dynamics evolve independently.Then, each eigenfunction corresponding to the zero eigenvalue is piecewise constant on one of the independent sectors.The optimal encoder Eq.D10 will depend instead on ϕ 0 (x), which identifies the independent sectors.

Appendix E: Rate of information decay
We consider the mutual information where x denotes values of the random variable X t and y the values of X t+∆t .We further consider dynamics which can be given by a transfer operator U with integral kernel p(y|x) that can be spectrally decomposed as H .The appearance of an unstable direction at β1 ≈ 3 coincides with the first IB transition.The emergence of a second unstable direction doesn't correspond precisely to the IB transition at β2 because stability is evaluated at the uniform encoder, but the true optimal encoder at β2 is not uniform.(c) The unstable directions are dominated by single components (note the color scale is logarithmic).At the first transition, the logarithm of the encoder is given by the eigenfunction ϕ1(x), up to rescaling (y-axis is shown in arbitrary units).(d) Likewise, at the second transition the encoder is given primarily by ϕ2(x).
where λ 0 = 0, ϕ 0 = const and ρ 0 (y) is the steady state distribution.The mutual information can then be written e (λn+λm)∆t x,y e (λn+λm+λ ℓ )∆t x,y For long times, the contribution of λ 1 dominates, with the remaining terms decaying as e (λn−λ1)∆t for n ≥ 2. Retaining only the first term, we see For short times, other eigenvalues will also contribute to the mutual information (Fig. S4).( 1 2 − ε) 2 terms subtracted with a constant (here 1/8).Here, the spectral gap λ0 − λ1 does not close but there is an accumulation of eigenvalues at ε → 0.5.(d-e) Mutual information for ∆t = 1 for the matrices in (b-c), respectively, as a function of ε.The curve from the other scenario is shown in gray for comparison.Both exhibit a clear peak in information which differs only slightly in magnitude, showing that the information peak in Fig. 2d-f of the main text is not due to the closing gap alone, but rather due to contributions from subdominant eigenvalues.and hence If p(h|x) and p(h) are chosen properly, the Kullback-Leibler can be expressed analytically which enables gradients to be effectively computed.As in Ref. [41], we take a Gaussian ansatz for p(h|x) and let the marginal p(h) be a spherical unit-variance Gaussian.More concretely, encoded variables H t are sampled from p(H t |X t ) by computing where f W and σ W are deterministic functions modeled by neural networks with parameters (weights) W , and η is a Gaussian random variable with unit variance.
To bound the entire loss from above, we must bound I(Y, H) from below.We do this using the noise-contrastive estimate of the mutual information introduced in Ref. [81].This recasts the problem as one of distinguishing samples from the distributions p(y|h) and p(y).Given a batch of B pairs (y, h) and one particular value h i , one asks what the probability is that a sample y j is from p(y|h i ) (is a positive sample) and not p(y) (is a negative sample).This probability is where i is the index of the positive sample.The log likelihood of these probabilities, where the expectation is taken over indices i and j, is closely related to the mutual information I(Y, H).In the limit of infinite samples B → ∞, this quantity is given by If one had access to the probabilities p(y i = pos|h) appearing in Eq.F7, the mutual information could thus be easily determined.One attempts to estimate these probabilities by introducing a variational ansatz f (y, h) to approximate the density ratio p(y|h) p(y) .Typically this f is represented by a neural network.One can then obtain a bound on the mutual information by minimizing from which the InfoNCE estimate of the mutual information can be calculated as The full objective to be minimized is given by an illustration can be seen in Fig. S5.This loss is evaluated on batches of sample pairs {(x t+∆t ), (x t+∆t ), ...}.A minimum is found via stochastic gradient descent.Note that other variational loss functions inspired by the information bottleneck have been derived, for example [82].
The VIB loss Eq.F11 is very similar to the loss function for a time-lagged β-variational autoencoder (β-VAE) [83], which have a form L = D KL − L rec , where the reconstruction term L rec measures the deviation from the true future state X t+∆t and the reconstruction from the latent variable H t .In the VIB, this term is replaced by the mutual information I NCE : rather than searching for a latent variable which can reconstruct the While a standard β-VAE with mean squared error (MSE) loss finds a latent variable that can reconstruct the full state at a time ∆t in the future, the VIB merely tries to reconstruct the statistics of X t+∆t conditioned on X t .
Replacing reconstruction losses with information-theoretical loss functions makes sense in some scenarios, such as chaotic systems.Recent work in this direction has used the Kullback-Leibler divergence as a loss function in place of a L 2 loss which generated well-behaved long-term dynamics [84,85].Other metrics which penalize the number of "false neighbors" in latent space have also shown to improve performance for chaotic systems [86].
The extent to which our findings might apply to time-lagged VAEs with L 2 reconstruction losses is an interesting question for future work.In Fig. S13, we observe that time-lagged VAEs learn a similar latent variable as VIB for the simulated fluids dataset.Indeed, several works have noted the apparent similarity between modes learned by variational and regular autoencoders with linear methods such as principal component analysis (PCA) [18,87,88].

VIB verification for noisy Hopf oscillator
We now show that VIB learns an encoding consistent with exact IB, i.e. one which depends only on the dominant transfer operator eigenfunctions.We consider the example of a Hopf oscillator given by the dynamical system Equation (F12) is the normal form (in polar coordinates) for dynamics near a Hopf bifurcation [89].For µ > 0 this system exhibits a circular limit cycle of radius √ µ (Fig. S5b).In Section G we show that this system has infinite purely imaginary eigenvalues λ n = inω and that as a result, encoders which exactly encode the angle coordinate, p(h|r, θ) ∝ δ(h − θ), are solutions of the IB optimization problem.The exact IB calculation breaks down for perfectly deterministic dynamics, hence we add a small amount of white noise to the dynamics.While this slightly perturbs the spectral content of the transfer operator, we still find our results to be consistent with the those expected for the deterministic case (see Section G).As we increase the size of the encoding alphabet N H , the encoder partitions space by finer and finer angular wedges.The learned encoding is invariant under permutations of the encoded symbols.Upon reordering, we find that, for large alphabets N H ≫ 1, the encoder indeed approximates p(h|r, θ) ∝ δ(h − θ) as expected (Fig. S5c).
Using VIB, we learn a continuous h t which can be computed directly from samples of the state variable.The encoding h t learned by VIB closely approximates the angle coordinate (Fig. S5e-f) and is nearly uncorrelated with r (inset).This shows that our mathematical results illustrated for exact IB in Sec.D hold also in the approximate framework of VIB.

VIB for deterministic dynamics
IB is an inherently probabilistic framework.To handle deterministic dynamics, we make them effectively stochastic by introducing a stochastic sampling scheme.Concretely, rather than taking as our IB relevance variable Y = X t+∆t , we take Y = X t+∆t+η where η is a random uniformly-distributed time shift.Despite using a different relevance variable, this yields essentially the same optimal encoding as (D1), where the eigenvalue e λn∆t is replaced by dη p(η) exp(λ n (∆t+η)).Crucially, the encoder retains its dependence on the transfer operator eigenfunctions ϕ n (x) as before.In general, this need not be the case: selecting a new relevance variable changes the IB objective and will generically lead to a different encoder.Our choice of stochasticity which we introduce to the dynamics is chosen in such a way to preserve the form of the encoder.
We illustrate this with a prototypical example of deterministic chaotic dynamics, the Lorenz system [90].In the steady state, the state variable of the Lorenz system resides on a chaotic attractor consisting of two "lobes" which encircle unstable fixed points (Fig. S6a).The encoding variable h learned by VIB is the slowest varying function of the state x, decaying as exp(λ Re 1 t) where λ Re 1 is the real part of the first subleading eigenvalue of the transfer operator (here Perron-Frobenius operator) (Fig. S6b).Eigenvalues λ n as well as right eigenvectors ϕ n (x) of the Perron-Frobenius operator were computed numerically using the Ulam method.The correspondence between ϕ 1 (x) and h disappears as β is increased (Fig. S6c-e).The dynamics of h also become notably less "slow", and are instead nearly constant with occasional large jumps (Fig. S6f).This shows that in order for h to be a valid slow variable, the compression term in Eq. (F11) is crucial.These results are consistent with those obtained for exact IB, where we showed that the encoder incorporates the slow modes at low beta (high compression).
Because this differential equation is separable, we can find eigenfunctions by looking for eigenfunctions that are a function of either r or θ.Recall that for deterministic dynamics, a product of adjoint transfer operator (or Koopman operator) eigenfunctions is again an eigenfunction [19].The eigenvalue equation for the radial coordinate is solved in Ref. [33], and in Ref. [76] it is shown that for this system there is no globally valid Koopman decomposition.However we are only interested in a subset of eigenfunctions, in particular those with eigenvalue with real part Reλ i ≈ 0, which may not suffice to approximate arbitrary functions of the state variable.
For the angle coordinate alone the situation is much simpler, which is solved by functions ϕ n (θ) = exp λn ω θ .Periodicity requires λ ω 2π = i2πn, which leads to λ = inω.The eigenfunctions and eigenvalues are then The corresponding eigenfunctions of the Perron-Frobenius operator are given by ρ = e −inθ .We now consider an encoding p(h|x).
where we retain only the eigenvalues with zero real part; the λ k in the above refer to those eigenvalues with non-zero real part.After neglecting these terms, it can be seen directly that p(h|r, θ) = p(h|θ).The coordinates (r ′ , θ ′ ) are used to denote (r t+∆t , θ t+∆t ).The expression can be further simplified by replacing the sum over n with a delta function, which leads to We next ask whether an encoder of the form p(h|θ) ∝ δ(h − θ) is a solution.Recall that the probability distribution appearing in Eqn.G9 is a distribution over future positions, p(R t+∆t = r ′ , Θ t+∆t = θ + ω∆t|h).This distribution can be calculated as where we used that the two terms in the integrand of the first equation are both delta functions.Plugging this into (G9) shows that p(h|θ) ∝ δ(h − θ) is consistent, and hence is an optimal encoding.In the presence of noise, the eigenvectors are perturbed and may gain a dependence on r (see Fig. S7).The presence of noise changes the eigenfunctions, in particular they depend on r.

Details on gradient analysis
The fluid flow of a von Karmen vortex street is well approximated by linear dynamics, which can be understood by recognizing that the system is poised just after a Hopf bifurcation so that there is an angle coordinate which rotates with constant angular velocity.This constant rotation can be described by a linear dynamical system.It follows from linearity that eigenfunctions of the adjoint transfer operator are given by linear functions of the state variable, where m (n) is the n-th "Koopman mode" and angled brackets denote integration over space.To compute these modes and the corresponding transfer operator spectrum, we use dynamic mode decomposition (DMD; see SI Section H) [11,12].We allow VIB to learn a two-dimensional encoding variable [h 0 , h 1 ], so that it can learn the complete first eigenfunction rather than only the real or imaginary part.For the purposes of comparing our learned variable with ϕ 1 , we construct a complex h = h 0 + ih 1 out of the two learned components.
To understand which function has been learned by the neural network, we check whether the learned functions are of the form Eq. (H1).This can be done by examining gradients of the network with respect to the input field where we have separated the part of the gradient which is independent of v from a residual part which is dependent on v. Gradients of the true eigenfunctions are given simply by We can then directly compare our latent variables with the true eigenfunctions by comparing their derivatives.If h corresponds to the true eigenfunction, we expect that m (IB) is approximately equal to the Koopman mode m (1) , and that g res is small.While it is unclear how to perform this decomposition in a general setting, we assume that the residual component g res,j averages to zero over an oscillation period of the flow field v.Then, m (IB) j ≈ ⟨ ∂h ∂vj ⟩ t and g res,j is given by variations about the mean.We see in Fig. S8 that these variations are much smaller than the mean in magnitude, and that they are essentially orthogonal to the mean vector.From this, we conclude that the gradients are given primarily by the constant part m.The learned latent functions vary for different training instances of the neural network.To extract the average gradient, we use PCA in an approach similar to that in Ref. [52].The learned functions can have arbitrary sign structure; h = h 0 + ih 1 is just as likely to be learned as h = −h 0 + ih 1 , for example.While in principle the network could learn arbitrary rotations, rather than simply changes in sign, we observe this is not the case.The distribution of gradients forms clusters in the high-dimensional gradient space corresponding to the four possible permutations of sign.PCA picks out the directions separating these clusters, as these are precisely the directions along which the data varies the most.This procedure gives an average gradient, while taking the varying sign structure into account.For reference, the gradient of a single instantiation of the VIB network can be seen in Fig. S13b.field, an order parameter can be computed as where V refers to the volume being integrated over.The value of r(t) is referred to as the synchronization order parameter, while ψ(t) is the average phase.To calculate the field θ(x, t) from an intensity field I(x, t) we imagine that the intensity field represents one component of the complex phase, for example I(x, t) = sin(θ(x, t)).The other component can be accessed using a time-delay, cos(θ(x, t)) = cos(θ(x, t + τ )) = I(x, t + τ ) =, where τ here should be chosen so that the intensity field undergoes one quarter of a full oscillation.As we know the true period of the circadian cycle is 24 hours/frames, we choose a delay of τ = 8 frames.Then, the phase can be computed as θ(x, t + τ ) ≈ arctan I(x, t) I(x, t + τ ) .
In Fig. S11 we show the results of VIB when applied to a locally-coupled Kuramoto model.Here we learn the same latent features which undergo oscillatory dynamics of varying radius, where the radius corresponds to the synchronization order parameter.

Appendix J: Simulation Parameters
Triple well simulations For the dynamics of the triple well we work directly with the force The evolution of the Brownian particle's position x t is given by where η t is white noise with unit variance.The noise magnitude σ is related to the diffusion constant in Eq.C1 by D = σ 2 2 .We use σ = 1.0 and γ = 0.2 Simulations were performed with N init = 10 5 initial conditions, with 300 trajectories generated from each initial condition.The state is evolved for 100 steps at dt = 2 • 10 −3 .The transfer matrix is approximated by binning the space x ∈ [−1, 1] with N bins = 100 bins.IB is performed using a time delay ∆t = 64 steps.Pitchfork bifurcation simulations The system was evolved according to the stochastic differential equation Eq.J1 with F (x) = −µx − x 3 , γ = 1 and σ = 0.1.For each value of µ, 10 5 initial conditions were simulated with 2000 trajectories starting at each initial condition.These were evolved for 100 time steps of dt = 2 • 10 −3 .Hopf oscillator simulations The system was evolved according to Eq. J1 with the force given by Eq.F12 in the main text, with γ = 1 and σ = 10 −2 , ω = 4.0 and µ = 0.25.For each value of µ, 10 6 initial conditions were simulated with 1000 trajectories starting at each initial condition.These were evolved for 50 time steps of dt = 2 • 10 −3 .Lorenz system simulations The system was evolved according to with ρ = 28, β = 8/3 and σ = 10.We took 10 3 initial conditions from which trajectories were simulated for 10 5 steps with dt = 2 • 10 −3 .We compute the true eigenfunctions using the GAIO library [92].
Fluid flow simulations The fluid flow simulations simulations are contained in the "Cylinder in Crossflow" dataset, downloaded from Ref. [93].We take the dataset at Reynolds number 150, and interpolate the velocity field from the unstructured 6569-node mesh to a regular grid of 300 × 150 pixels.The VIB networks are trained with β = 10 7 and a time buffer ∆t chosen randomly between [4,24] (see our comment on time randomization in "VIB for deterministic dynamics" in Section F).
Locally-coupled Kuramoto model In SI Fig. S11 we perform variational IB on a locally-coupled Kuramoto model and find the latent variables learn a synchronization order parameter as in the cyanobacteria data.This model considers a spatially-distributed field of phases which evolves according to where x i denotes the location of gridpoint i and N i denotes the sites neighboring i.These dynamics are integrated using the Euler method, where we take ω(x) = const = 0.05 for the natural frequencies, and J = 1.0.

Appendix K: Comparison of various dimensionality reduction methods
In this section we compare the performance of variational IB (VIB) to several common data-driven model reduction or inference methods.VIB uses a neural network to identify the relevant variables that are most predictive of the future.This construction learns dynamically relevant variables in contrast to methods such as principal component analysis or diffusion maps (in the absence of additional stringent assumptions), and it finds potentially non-linear variables in an agnostic way without requiring a tailored library of non-linearities (such as extended dynamic mode decomposition).In addition to these factors, VIB produces physically well-defined relevant variables -transfer operator eigenfunctions -in contrast to other deep methods such as (variational) autoencoders.In Fig. S12 we compare the low-dimensional latent space trajectories of the reduced models obtained using several of these standard methods to those output by VIB.In the remainder of this section we discuss these methods in detail.

Principal Component Analysis/Proper Orthogonal Decomposition
Principal component analysis, or PCA, is a linear method that projects the data onto a subspace which accounts for the most variance in the dataset [94].In the dynamical systems literature, PCA is also known as a proper orthogonal decomposition (POD) [95].The starting point for PCA is a collection of samples x (k) , where each sample has N feat different features: x (k) ∈ R N feat .Given these samples, one computes a correlation matrix between every pair of features The principal components of the data are eigenvectors of the covariance matrix; the dominant principal component determines the direction of most variation in the dataset.In practice, if many samples are present the correlation matrix is expensive to compute.Writing the dataset as a matrix, X ∈ R N samples ×N feat , computing the correlation matrix requires evaluating the matrix product C = X T X.Instead, a more computationally efficient way is to compute the singular value decomposition (SVD) of the data matrix, X = UΣV T .
The principal components are given by columns of V, and the singular values in the matrix Σ are the square root of the variance along each principal component direction.PCA's relation to SVD also means that one can interpret a reconstruction of X using only k principal components as the optimal rank-k approximation of the dataset, where optimality is here defined with respect to the Frobenius norm.
Compared to VIB, the linearity assumptions underlying PCA can be viewed either as a benefit or a drawback.The benefit is that a linear projection can be intuitive to understand in terms of the original data, and that one can compute this projection quickly.The potential downside is that a linear projection may not be sufficient or appropriate for systems that are highly non-linear.In addition, PCA makes no reference to the dynamics of the system, unlike VIB.An example of where PCA can fail is shown in Fig. 3 of Ref. [73].In words, imagine a particle is hopping in a 2D double-well potential, where the wells are centered at x ± 1.If the vertical extent of the wells is very large (|y| ≫ 0), then the dominant principal component will be along the y-direction, even though the interesting dynamics are the hops between wells, which happens in the x-direction.

Dynamic Mode Decomposition
Dynamic Mode Decomposition (DMD) starts from the assumption that the system evolves linearly [11,12].Concretely, if x t ∈ R n feat denotes the observed quantity at one point in time, then one assumes dynamics given by The matrix A is a finite-dimensional approximation to the Koopman operator [19].
To find the matrix A DMD , we start by assembling a collection of samples into a data matrix X t ∈ R N samples ×N feat .In addition to this, we assemble a time-shifted matrix X t+1 of the same shape as X t .In the scenario where we only have one trajectory of duration T , each sample may correspond to the system at a measured time point (excluding the final time), so that The matrix A DMD is given by the least squares solution where M + denotes the pseudoinverse of the matrix M. Eigenvectors w of the matrix A DMD are Koopman modes, which correspond to left eigenfunctions of the transfer operator.The evolution of these Koopman modes in time is given by the time-dependent amplitude a i (t) = x t • w, where we take x t here to be a single measurement at time t.
Similar to PCA, DMD reduces the dimensionality of the dataset by finding a linear projection of the data onto a subspace of the full features space.However, a key difference between the two approaches is that DMD uses the system's dynamics to identify an "optimal" subspace, whereas PCA identifies a subspace based only on the steady-state distribution of the data.
DMD in its original formulation assumes that the Koopman operator linearly evolves the observed state variable.However, the true Koopman operator evolves arbitrary non-linear functions of the state variable forward in time.DMD has therefore been extended to account for non-linear functions through "extended DMD", or eDMD [13].eDMD augments the state vector with non-linear transformations of the state.As an example, a state vector [−x−] may be replaced by monomials [−x−, −x 2 −, −x 3 −] (where here we understand exponentiation as element-wise).Other approaches also exist, such as augmenting the state vector by several time-delayed state vectors [96].Rather than constructing the full d × d matrix, where d is the dimension of the (possibly augmented) state, one can directly construct a low-rank approximation of K by computing the reduced-rank singular value decomposition (SVD) of the matrix K from the SVDs of X t and X t+∆t .
Similar to PCA, DMD is primarily limited in its assumption of linear dynamics.In some cases this can be resolved with eDMD, which requires that one identifies a suitable set of non-linear terms to account for potential non-linear eigenfunctions of the Koopman operator.It is unclear how to choose this set of functions in a generic setting, which constitutes the biggest disadvantage to VIB.The features for eDMD must be hand-selected, which in some cases may defeat the purpose of using it as a feature-learning tool.VIB is not subject to this restriction, and can learn relevant features directly from the data.

Independent Component Analysis
Independent component analysis (ICA) was originally formulated in the context of blind source separation, where a useful picture is the "cocktail party problem" [10].Imagine you place a set of microphones in a room at a cocktail party; these microphones will record the combination of all the conversations happening at once.The goal of blind source separation is to find a way to, from the recorded signals, isolate the original conversations.Mathematically, ICA finds a solution to the equation Here, x denotes the recorded signal, s denotes the independent sources (conversations), and A ICA is the mixing matrix.The only measured quantity is the vector x; both the mixing matrix and the independent sources must be learned.To find one independent component, we start with an initial guess y = b T x which, for a correct choice of b, should be equal to some component s i .To identify the correct b we will use the fact that sums of independent variables are more Gaussian than the original variables themselves, which follows from the central limit theorem.By assumption, the observed signals are a linear combination of independent components, so y is also: y = q T s.If multiple q i are non-zero, this will be more non-Gaussian than if only one is non-zero.Thus, y = s i for the choice of b for which y is maximally non-Gaussian [10].
Numerically, one optimizes the deviation from Gaussianity using an approximation of the "negentropy" of the distribution of y: J neg = S(y Gaussian ) − S(y), where y Gaussian is a normally distribution random variable and S is the Shannon entropy.The intuition for this formulation is that the J neg is minimized if y has a unit-Gaussian distribution, so that it serves as a metric for non-Gaussianity.In practice one uses an approximation for J neg , see Ref. [10] for details.The optimal b can then be found by doing gradient ascent on this objective function.
On some level, independent component analysis (ICA) is similar to PCA in that it searches for a linear projection of the data onto a subspace of the feature space.In its basic implementation ICA, like PCA, does not incorporate dynamics in contrast to DMD or VIB.The essential difference between PCA and ICA is the assumption of independence of the components.In PCA, unless the data distribution is actually multivariate Gaussian, the components are unlikely to be independent.IB makes no assumptions about the independence of learned encoding variables with each other, only about their level of independence with the original state.We note that ICA can be understood as minimizing the mutual information between encoding variables, I(H 1 ; H 2 ), compared to IB's minimization of the mutual information with the original state [10].Which objective is more desirable depends on the system at hand, but in general we do not expect these to lead to the same encodings.

Time-lagged independent component analysis
The leveraging of non-Gaussianity in ICA is required due to the fact that, if using whitened data (z = Vx so that z is zero mean ,E[z] = 0, and unit variance, E[zz T ] = I), then the mixing matrix A ICA cannot be estimated if the components s are normal Gaussian variables: any orthogonal matrix will satisfy (K1).For dynamical systems, however, one can get around this requirement of non-Gaussianity.In particular, one can use a time-lagged variable z t−τ = A ICA s t−τ to compute the mixing matrix from the time-correlation matrix of the signal z t [10].In particular, we have where D is a diagonal matrix.In going to the bottom line, we used the assumption that the components s i are independent not just instantaneously, but also for a time lag τ .From here we can read off that the matrix A ICA is composed of eigenvectors of the correlation matrix of the whitened signal.When written in terms of the original (unwhitened) data x t and x t−τ it can be shown these eigenvectors are nothing other than the eigenvectors of the transfer matrix A DMD [73].Thus, the two methods are equivalent.

Diffusion maps
Diffusion maps are a technique which attempts to approximate the Perron-Frobenius operator, or rather the integral kernel p(x t+∆t |x t ) [14,97].Given this approximation, one computes eigenfunctions and uses them as a low-dimensional parameterization of the data ("diffusion coordinates").
This method takes data pairs {x t+∆t } i and approximates the probability of observing these two points via a kernel p(x t+∆t ), where one typically takes a Gaussian Ansatz From this, one can assemble the conditional probability distributions into a matrix .
This matrix describes the evolution of probability distributions on a graph where each node is a data point.In practice, a symmetrized version of P is constructed and the learned eigenvectors are adjusted after the diagonalization [14,97].
To compute the diffusion coordinates of an arbitrary point that wasn't in the original dataset, one inverts the definition of the adjoint transfer operator eigenfunction where λ i is the i-th eigenvalue.Diffusion maps have the advantage, relative to DMD, that they find the full Perron-Frobenius operator and not a linear approximation to it.However, while DMD can isolate the dominant eigenvectors of the operator using reduced-rank SVD, it is less clear how they can be extracted with diffusion maps without first computing the full matrix P .We note that VIB, like DMD, also directly learns the dominant modes and does not require estimation of the full transfer operator.

Deep (Variational) Autoencoders
Autoencoders (AEs) belong to a class of deep learning methods used for model reduction.Such approaches have successfully been applied in various domains for forecasting the dynamics of complex systems in terms of simpler latent dynamics [98,99].Autoencoders are composed of an encoder which compresses the observable x into a lower-dimensional latent variable z, and a decoder which attempts to reconstruct the original state x.Variational autoencoders aim to learn a probability distribution over observations p θ (x) ≈ p(x) from which one can directly sample [100].This is done by assuming that the latent variable is low dimensional, and optimizing the objective L β-VAE = E q ϕ (z|x) [− log p θ (x|z)] − βD KL (q(z|x)∥p(z)) where q ϕ (z|x) is the posterior on the latent variables z and is parameterized by a neural network (encoder) with parameters ϕ, and p θ (x|z) denotes the decoding network with parameters θ.Strictly speaking, we present in the above objective function the β-VAE loss [83].With β = 1, which is the case for the original VAE [100], the objective is an upper bound on the log likelihood −E[log p θ (x)] which is minimized when p θ (x) is equal to the true data distribution p(x).The term β controls compression as in the IB objective, however it has the opposite effect: small β IB corresponds to high compression, while small β VAE corresponds to low compression.
In cases where one directly computes the probabilities p θ (x|z) the first term can be evaluated directly, else it is typically replaced with an L 2 loss ∥x − g θ (z)∥ 2 (where g θ is a deterministic neural network), which is equivalent to assuming a Gaussian Ansatz for p θ with a fixed variance.
The VIB loss function is very similar to the β-VAE loss [101].Rather than attempting to reconstruct the original state x, VIB replaces this term with an estimate of the mutual information between the latent variable z and some other relevance variable y (for us, y = x t+∆t ).In contrast to DMD and diffusion maps, neither VIB nor β-VAEs make any mention of transfer operators and are instead motivated by purely statistical considerations.As we show in the main text, the latent variables learned by VIB correspond to eigenfunctions of the transfer operator.While we cannot claim that that β-VAEs learn the same thing, some preliminary results in Fig. S13 suggest they may coincide to some degree.In the cyanobacteria dataset, we see that the latent variables learned by a β-VAE show the same qualitative structure as in VIB, but are less smooth (Fig. S12).

Other neural networks
Other deep architectures can also be used for model reduction.For example, Ref. [23] uses recurrent neural networks (RNNs) to learn the evolution of macroscopic variables.To ensure stability and fidelity, the macroscopic variables are periodically "lifted" to the full microscopic state, which is then evolved for several time steps to recalibrate the RNN's hidden state.While interpretability of the latent variables has yet to be explored in such models, we expect the addition of VIB-like objective functions may aid interpretability without harming performance.
As another example, [20] attempts to learn the full Koopman operator using neural networks.This can be thought of as an extension of eDMD, where instead of prescribing the library of nonlinear terms by hand, they can be learned by a neural network.The latent dynamics are then encouraged to be linear by minimizing the deviation between the true future (latent) state z t+∆t and its linear approximation found by least squares.Similar approaches have also been explored in Refs.[21,22].As an illustration, we trained one such network on the cyanobacteria which can be seen in Fig. S12.This approach was applied primarily to deterministic systems.We expect that combining their method of encouraging linear latent dynamics together with the VIB objective function may be fruitful and lead to more well-behaved latent variables.

FIG. 2 .
FIG. 2. Information loss of a Brownian particle in a double well potential.(a) Fixed point (FP) diagram of the dynamics given by Eq. 7 for zero noise.There is a bifurcation at µ = 0 where the stable FP at x = 0 becomes unstable and two new stable FPs appear at ± √ µ.Insets show the evolution of the corresponding potential V (x), with the emergence of a double-well structure for µ > 0. (b) Dynamics of the system Eq.7 for varying values of µ corresponding to the potential insets in (a), with uniformly-distributed initial conditions.(c) Loss of information between the initial condition and the future state.Inset shows scaling given by the first eigenvalue of the transfer operator.(d) Mutual information between the present and future state for varying time delay ∆t and bifurcation parameter µ.(e) Spectrum of the transfer operator U , showing a pile-up of eigenvalues for µ ≳ 0. These are related to the eigenvalues of its infinitesimal generator by Λi = e λ i ∆t .(f) Maximal mutual information which can be encoded into a discrete variable of NH values, for a fixed time delay ∆t = 1.0.Black dashed line shows I(Xt, Xt+∆t) for reference.Information is provided in units of bits.

FIG. 3 .
FIG.3."Knowing when to stop".The spectral properties of the transfer operator determine the necessary complexity (i.e."when to stop"[25]) of the reduced model, which we show is also visible in information theoretic metrics.(a) Four-well potential in which a Brownian particle fluctuates.The magnitude σ of the fluctuating noise is related to an energy scale Eσ = σ 2 .(b) Information contained in the encoding variable Ht about the future state Xt+∆t for varying levels of noise and alphabet sizes NH .(c) Information gain achieved by increasing the alphabet size by a single variable.This is the discrete derivative of the curve in (b).(d) Spectrum of the transfer operator for changing values of noise amplitude.

FirstFIG. 4 .
FIG. 4. IB learns eigenfunctions of the adjoint transfer operator.(a) When the relative weight β between both constraints in Eq. 1 is changed, more and more information can go through the encoder.This occurs in steps, where the spectral content of the transfer operator is included starting from eigenvalues with largest magnitude (i.e., the slowest ones).(b) Transitions are characterized by the appearance of negative eigenvalues in the spectrum of the Hessian of the IB loss function.Here we consider the Hessian evaluated at the uniform encoder p(h|x) = N −1 H .The IB transitions β1 ≈ 1.1 and β2 ≈ 1.8 correspond to the appearence of negative eigenvalues of the Hessian.(c) The unstable directions are dominated by single components (note the color scale is logarithmic).(d) At the first transition, the logarithm of the encoder is given by the eigenfunction ϕ1(x), up to rescaling (y-axis is shown in arbitrary units).(e) Likewise, at the second transition the encoder is given by ϕ2(x).

FIG. 5 .
FIG. 5. Variational IB for high-dimensional simulated fluid flow.(a) A fluid flows into the system with uniform velocity v0 in the x-direction and passes a disk-shaped obstacle, which perturbs the fluid and causes vortex shedding behind the object in a so-called von Kármán street.The state of the system is given by a spatially varying two-component vector field v(x).(b) The dynamics in latent space (blue) are very regular, traversing a nearly circular trajectory (see Supplementary Movie 1).For comparison we show the evolution of the mode amplitudes obtained by projecting the velocity field onto the first DMD mode (green).(c) Time evolution of one component of the latent variable (h1, blue) as well as the DMD mode amplitude (green).(d) Comparison of the first Koopman mode obtained from DMD (m(1) ) and from VIB (m (IB) ).Koopman modes from VIB are computed as gradients of the latent encoding variables as described in the main text.Red corresponds to positive values and blue to negative; the magnitudes of the modes are not directly comparable.

FIG. 6 .
FIG. 6. Discovering slow collective variables in cyanobacteria populations (a) Fluorescent images of cyanobacteria colonies labeled with EYFP driven by the kaiBC promoter, allowing the visualization of Kai protein transcription.The colonies are imaged as they undergo cell growth and oscillations in Kai expression associated with the circadian rhythm.(b) Time series of one cropped section of an individual colony (see SI Fig. S9 for details).(c) The "state" used for variational IB is the time-lagged intensity field with lag time τ .(d) Variational IB embeddings of time-lagged images into two dimensions (see Supplementary Movie 2).Every line corresponds to one colony's evolution in time.Note the apparent oscillations of different radii.(e) Embedding into three dimensions.We orient our axes to correspond to the three principal components of the data.Thus, the projection onto the h2 = 0 plane corresponds to a projection along the dominant two principal components.In this subspace, we see a similar structure as in the 2D embedding.(f) Mutual information between the future state and the encoding, given by I(Xt+∆t, Ht), for varying dimension of the latent space.Increasing the dimension of the embedding space beyond two leads only to marginal increases in I(Xt+∆t, Ht); this tells us "when to stop."Each small point represents one training instance of the variational IB model, while the large point shows the maximum estimated value.Because the InfoNCE estimator is a lower bound on the true information, we consider the maximum as the estimated mutual information.The value of Imax = I(Xt, Xt+∆t) is the mutual information estimated for the true dynamics.(g) Selected trajectories in latent space with large and small radii.(h) Time-average latent radius of all colonies.(i) Mean synchronization order parameter of the intensity images; see SI for computation details. (j) Mean radius versus synchronization parameter for each colony.VIB identifies clusters of cells characterized by high (blue) or low (black) synchronization.As revealed to us after our analysis, this clustering corresponds to differing theophylline concentrations across experiments.Within each experimental movie (each of which contains 2-3 colonies) the radii are mostly constant.(k) Sample time series of weakly (black) and highly (blue) synchronized colonies.We apply a slight Gaussian blur to better visualize the bacteria boundaries.
FIG. S2.Information bottleneck for a Brownian particle in a double well.(a) Mutual information between the current state Xt and its encoding Ht for various alphabet sizes NH (color).The first two IB transitions are denoted with black lines and occur at β1 and β2.(b) Mutual information between the current state's encoding Ht and the future Xt+∆t.The gray dashed line denotes the maximum attainable information, which is the mutual information I(Xt, Xt+∆t).Color denotes encoding alphabet size.Black lines showing IB transitions are shown for reference.(c) Optimal encoder after the first transition (β ≳ β1) with alphabet size NH = 3.(d) Optimal encoder after the second transition (β ≳ β2) with alphabet size NH = 3.
FIG. S3.IB transitions for a Brownian particle in a triple-well potential.(a) Information transitions with varying β for a particle in a triple-well potential.(b) Eigenvalues of the IB Hessian evaluated at the uniform encoder p(h|x) = N −1H .The appearance of an unstable direction at β1 ≈ 3 coincides with the first IB transition.The emergence of a second unstable direction doesn't correspond precisely to the IB transition at β2 because stability is evaluated at the uniform encoder, but the true optimal encoder at β2 is not uniform.(c) The unstable directions are dominated by single components (note the color scale is logarithmic).At the first transition, the logarithm of the encoder is given by the eigenfunction ϕ1(x), up to rescaling (y-axis is shown in arbitrary units).(d) Likewise, at the second transition the encoder is given primarily by ϕ2(x).
FIG. S4.Information gain due to an eigenvalue pile-up.(a) To study the role of a gap closing compared to the pile-up of eigenvalues (beyond the dominant one) in the mutual information I(Xt; Xt+∆t) we build a matrix P with which we can tune the eigenvalues via a parameter ε.(b) Eigenvalues for the matrix shown in (a) with varying parameter values ε; for ε → 0.5, the spectral gap λ0 − λ1 closes.(c) Eigenvalues for the matrix shown in (a), but with all 12 ( 1 2 − ε) 2 terms subtracted with a constant (here 1/8).Here, the spectral gap λ0 − λ1 does not close but there is an accumulation of eigenvalues at ε → 0.5.(d-e) Mutual information for ∆t = 1 for the matrices in (b-c), respectively, as a function of ε.The curve from the other scenario is shown in gray for comparison.Both exhibit a clear peak in information which differs only slightly in magnitude, showing that the information peak in Fig.2d-f of the main text is not due to the closing gap alone, but rather due to contributions from subdominant eigenvalues.
FIG. S5.Variational IB extends exact IB.(a) Variational IB optimizes the variational objective directly from samples, in contrast to exact IB which requires an estimation of the full conditional distribution p(xx+∆t|xt).Values of ht are drawn from a "latent" distribution p(h|x), from which one can estimate the mutual information I(Xt+∆t, Ht).The compression term is approximated by the Kullback-Leibler divergence between the learned p(h|x) and a variational ansatz for the marginal p(h).(b) Phase portrait for a dynamical system above a Hopf bifurcation.(c) Exact IB learns (up to permutations) an encoding h which corresponds to the polar angle coordinate θ.(d) Variational IB similarly learns the correct encoding, directly from samples.(e) Correspondence between IB encodings and the angle coordinate θ.The encoding learned by VIB is independent of r (inset).(f) Time dependence of the encoded variable ht.

FIG. S6 .
FIG. S6.Variational IB discovers slow variables.(a) Variational IB applied to the Lorenz system.Color corresponds the encoding value ht.(b) Time correlation functions of the encoded variable CXY (τ ) = ⟨(X(τ )− X)(Y (0)− Ȳ )⟩ for high compression (black) and no compression (red).Blue shows the decay of the subleading Koopman eigenfunction using numerically-obtained eigenvalues.(c-d) Lorenz attractor projected onto the x-y plane, colored by encoding variable.(e) Same projection, colored by the value of the true subleading Koopman eigenfunction.(f) Dynamics of ht obtained by encoding the state at every time during the true trajectory.Top corresponds to no compression (β = ∞), bottom corresponds to high compression (β = 4000).
FIG. S7.Eigenfunctions and IB partitions for the Hopf Oscillator.(a) Phase portrait for the deterministic Hopf oscillator.(b) Partitions found by IB for high beta but a restricted encoding alphabet NH .(c) Example trajectories of the nearly deterministic Hopf oscillator, with small noise amplitude σ = 0.01.(d) IB partition of the low noise dynamics for a NH = 16 encoding alphabet (same as appears in (b)).(e) For small noise, the first several eigenfunctions obtained numerically approximate the numerically expected ones of the form cos nθ.(f) Simulated trajectories of Hopf oscillator with higher noise amplitude σ = 0.2.(g) IB partition of the low noise dynamics for a NH = 16 encoding alphabet; note the dependence on r.(h) The presence of noise changes the eigenfunctions, in particular they depend on r.
FIG.S8.VIB gradients for fluid flow are nearly linear.We compare the linear part of the gradients of the VIB network to the residual which depends on the input field v(x).(a) The residual parts are smaller in magnitude than the linear part m.Angle brackets denote average over space.(b) The residual parts are nearly orthogonal to the mean m, as we see that the projection of the full gradient onto the mean is nearly 1.
FIG. S9.Cyanobacteria data and choice of time delay (a)The cyanobacteria dataset consists of 5 videos, where each has multiple (up to 3) colonies.(Top) time evolution of one video, where the interiors of the colonies are outlined with white boxes.The dataset is composed of all colonies which are greater than 64px in height and width.(Bottom) evolution in time of one selected colony.(b) Learned VIB encodings when no time delay is used, i.e. τ = 0. Left shows the evolution of (a subset of) individual colonies, colored by the original experimental video which they belong to.Note that the axes scale is adjusted for each trajectory so that it fills the plot.Right panel shows all colonies in latent space.(c) Same as above, but with the 3 hr time delay used in the main text.(d) Hyperparameter sweep over different choices of time delay τ and prediction horizon ∆t.Here each set of three plots was trained with the same parameters but different neural network instantiations to understand how robust these latent spaces are.
FIG. S10.Variance along principal component directions For each DVIB model trained on the cyanobacteria dataset, we compute the principal component decomposition of the resulting point cloud in latent space.Here we show the variance along each principal component for 60 instantiations of the DVIB model.Variances are normalized by the variance in the first principal component direction.
FIG. S11.VIB applied to simulated Kuramoto model Inputs are made to mimic the intensity field of the cyanobacteria data.The intensity field (a) is derived from the phase field (b) according to I(x, t) = 1 2 cos θ(x, t) − 1. (c) DVIB trained with two latent variables learns oscillations.Here each trajectory is colored by its synchronization order parameter.(d) Latent oscillations are very highly correlated with synchronization order parameter.
FIG. S12.Comparison of various dimensionality reduction techniques on cyanobacteria dynamicsTop row illustrates linear methods, bottom row shows non-linear methods.For each method, we plot the trajectories of all cyanobacteria colonies when projected onto the dominant modes (whose definition depends on the particular method).Blue trajectories show the highly synchronized colonies corresponding to high theophylline concentration, black/gray trajectories show the unsynchronized colonies from low theophylline conditions.Unlike the first four methods, deep neural networks (bottom right) directly produce a two-dimensional latent variable and do not provide access to subleading modes.(a) Principal component analysis reduces the dimensionality of data by finding a projection onto the directions of most variance (sketch).Here we show the projection of the cyanobacteria trajectories onto several pairs of principal components, ai(t) = x(t) • mi.Each curve corresponds to the trajectory of one single colony evolving under one of four theophylline conditions; there are 10 trajectories in total, coming from 4 experimental conditions.Below, the first four principal components are shown.Here, the state is composed of two time-lagged images of the cyanobacteria colony; the top half of m corresponds to the earlier image, the bottom half corresponds to the lagged image.(b) Dynamic mode decomposition is performed on pairs of data points and aims to find a linear evolution operator A that links them (sketch).(Right) Trajectories of the projection of the state onto several eigenvectors of A. The eigenvectors themselves are shown below (the imaginary parts of the first two modes are uniformly zero).(c) Independent component analysis seeks a linear projection of the data onto statistically independent components si.In other words, the distribution p(si|sj) is independent of sj (sketch).(Right, top) Trajectories of the projection of the state onto several independent components.(Bottom) columns of the mixing matrix; the components si describe how the weighting of these columns evolves in time.(d) (Left) Diffusion maps build a graph of the observed state variables with edges weights determined by their distance from one another.(Right) Projection of the trajectories onto the first non-trivial eigenfunctions (because a0 =const) of the inferred transfer operator on the graph.(e) A (time-lagged) autoencoder is a deep neural network architecture which encodes the state xt into a low-dimensional latent space.From the encoding one then tries to reconstruct the future state.(Right) Trajectories of the cyanobacterial colonies when encoded into a two-dimensional latent space.(f) A variational autoencoder is similar to a standard autoencoder, where one instead learns a distribution over possible encodings, p(ht|xt).(Right) Trajectories of the cyanobacterial colonies when encoded into a two-dimensional latent space.(g) Implementation of the network in Ref.[20] which encodes the state using a time-lagged autoencoder with an additional loss term that aims to enforce linearity of the embedded dynamics, ht+∆t ≈ Kht.(h) Latent trajectories produced by variational IB (VIB).