Foundation Inference Models for Markov Jump Processes

Authors: david  berghaus
kostadin  cvejoski
patrick  seifner
césar  ojeda
ramsés j sánchez
lamarr  institute
fraunhofer  iais
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-13 17:03:14

Abstract: Markov jump processes are continuous-time stochastic processes which describe dynamical systems evolving in discrete state spaces. These processes find wide application in the natural sciences and machine learning, but their inference is known to be far from trivial. In this work we introduce a methodology for *zero-shot inference* of Markov jump processes (MJPs), on bounded state spaces, from noisy and sparse observations, which consists of two components. First, a broad probability distribution over families of MJPs, as well as over possible observation times and noise mechanisms, with which we simulate a synthetic dataset of hidden MJPs and their noisy observations. Second, a neural recognition model that processes subsets of the simulated observations, and that is trained to output the initial condition and rate matrix of the target MJP in a supervised way. We empirically demonstrate that *one and the same* (pretrained) recognition model can infer, *in a zero-shot fashion*, hidden MJPs evolving in state spaces of different dimensionalities. Specifically, we infer MJPs which describe (i) discrete flashing ratchet systems, which are a type of Brownian motors, and the conformational dynamics in (ii) molecular simulations, (iii) experimental ion channel data and (iv) simple protein folding models. What is more, we show that our model performs on par with state-of-the-art models which are trained on the target datasets.

Our pretrained model is available online.

Main paper: 
 introduction: 
 Very often one encounters dynamic phenomena of wildly different nature, that display features which can be reasonably described in terms of a macroscopic variable that jumps among a finite set of long-lived, metastable discrete states. Think, for example, of the changes in economic activity of a country, which exhibit jumps between recession and expansion states (Hamilton, 1989) , or the internal motion in proteins or enzymes, which feature jumps between different conformational states (Elber and Karplus, 1987) . The states in these phenomena are said to be long-lived, inasmuch as every jump event among them is rare, at least as compared to every other event (or subprocess, or fluctuation) that composes the phenomenon and that occurs, by construction, within the metastable states. Such a description in terms of macroscopic variables effectively decouples the fast, intra-state events from the slow, inter-state ones, and allows for a simple probabilistic treatment of the jumping sequences as Markov stochastic processes: the Markov Jump Processes (MJPs). In this work we are interested in the general problem of inferring the MJPs that best describe empirical (time series) data, recorded from dynamic phenomena of very different kinds. To set the stage, let us assume that we want to study some D-dimensional empirical process z(t) : R + → R D , which features long-lived dynamic modes, trapped in some discrete set of metastable states. Let us call this set X . Let us also assume that we can obtain a macroscopic, coarse-grained representation from z(t) -say, with a clustering algorithm -in which the fast, intra-state events have been integrated out (i.e. marginalized). Let us call this macroscopic variable X(t) : R + → X . If we now make the Markov assumption and define the quantity f (x|x ′ )∆t as the infinitesimal probability of observing one jump from state x ′ (at some time t), into a different state x (at time t + ∆t), we can immediately write down, following standard arguments (Gardiner, 2009) , a differential equation that describes the probability distribution p MJP (x, t), over the discrete set of metastable states X , which encapsulates the state of the process X(t) as time evolves, that isEQUATION)Equation 1 is the so-called master equation of the MJP whose solutions are completely characterized by an initial condition p MJP (x, t = 0) and the transition rates f : X × X → R + .With these preliminaries in mind, we shall say that to infer an MJP from a set of (noisy) observations z(τ 1 ), . . . , z(τ l ) on the process z(t), recorded at some observation times τ 1 , . . . , τ l , means to infer both the transition rates and the initial condition determining the hidden MJP X(t) that best explains the observations. In practice, statisticians typically assume that they directly observe the coarsegrained process X(t). That is, they assume they have access to the (possibly noisy) values x 1 , . . . , x l , taken by X(t) at the observation times τ 1 , . . . , τ l (see Section 2). We shall start from the same assumptions. Statisticians then tackle the inference problem by (i) defining some (typically complex) model that encodes, in one way or the other, equation 1 above; (ii) parameterizing the model with some trainable parameter set θ; and (iii) updating θ to fit the empirical dataset.One issue with this approach is that it turns the inference of hidden MJPs into an instance of an unsupervised learning problem, which, as history shows, is far from trivial (see Section 2). Another major issue is that, if one happens to succeed in training said model, the trained parameter set θ * will usually be overly specific to the training set {(x 1 , τ 1 ), . . . , (x l , τ l )}, which means it will likely struggle to handle a second empirical process, even if the latter can be described by a similar MJP. Figure 1 contains snapshots from two empirical processes of very different nature. The figure on the left shows a set of observations (blue circles) recorded from the discrete flashing ratchet process (black line). The figure on the right shows the ion flow across a cell membrane, which jumps between different activity levels (blue line). Despite the vast differences between the physical mechanisms underlying each of these processes, the coarse-grained representations of the second one (black line) is abstract enough to be strikingly similar to the first one. Now, we expect that -at this level of representation -one could train a single inference model to fit each process (separately).Unfortunately, we also expect that an inference model trained to fit only one of these (coarse-grained) processes, will have a hard time describing the second one.In this paper we will argue that the notion of an MJP description (in coarse-grained space) is simple enough, that it can be encoded into the weights of a single neural network model. Indeed, instead of training, in an unsupervised manner, a complex model (which somehow encodes the master equation) on a single empirical process; we will train, in a supervised manner, a simple neural network model on a synthetic dataset that is composed of many different MJPs, and hence implicitly encodes the master equation. This procedure can be understood as an amortization of the probabilistic inference process through a single recognition model, and is therefore akin to the works of Stuhlmüller et al. (2013) , Heess et al. (2013) and Paige and Wood (2016) . Rather than treating, as these previous works do, our (pretrained) recognition model as auxiliary to Monte Carlo or expectation propagation methods, we employ it to directly infer hidden MJPs from various synthetic, simulation and experimental datasets, without any parameter fine-tuning. We thus adopt the "zero-shot" terminology introduced by Larochelle et al. (2008) , by which we mean that our procedure aims to recognize objects (i.e. MJPs) whose instances (i.e. noisy and sparse series of observations on them) may have not been seen during training. We have recently shown that such an amortization can be used to train a recognition model to perform zero-shot imputation of time series data (Seifner et al., 2024) . Below we demonstrate that it can also be used to train a model of minimal inductive biases, to perform zero-shot inference of hidden MJPs from empirical processes of very different kinds, which take values in state spaces of different sizes. We shall call this recognition model Foundation Inference Model2 (FIM) for Markov jump processes.In what follows, we first review both classical and recent solutions to the MJP inference problem in Section 2. We then introduce the FIM methodology in Section 3, which consists of a synthetic data generation model and a neural recognition model. In Section 4 we empirically demonstrate that our methodology is able to infer MJPs from a discrete flashing ratchet process, as well as from molecular dynamics simulations and experimental ion channel data, all in a zero-shot fashion, while performing on par with state-of-the-art models which are trained on the target datasets. Finally, Section 5 closes the paper with some concluding remarks about future work, while Section 6 comments on the main limitations of our methodology. 
related work: 
 The inference of MJP from noisy and sparse observations (in coarse-grained space) is by now a classical problem in machine learning. There are three main lines of research. The first (and earliest) one attempts to directly optimize the MJP transition rates, to maximize the likelihood of the discretely observed MJP via expectation maximization (Asmussen et al., 1996; Bladt and Sørensen, 2005; Metzner et al., 2007) . Thus, these works encode the MJP inductive bias directly into their architecture.The second line of research leverages a Bayesian framework to infer the posterior distribution over the transition rates, through various Markov chain Monte Carlo (MCMC) algorithms (Boys et al., 2008; Fearnhead and Sherlock, 2006; Rao and Teg, 2013; Hajiaghayi et al., 2014) . Accordingly, these simulation-based approaches encode the MJP inductive bias directly into their trainable sampling distributions. The third one, also Bayesian in character, involves variational inference. Within it, one finds again MCMC (Zhang et al., 2017) , as well as expectation maximization (Opper and Sanguinetti, 2007) and moment-based (Wildner and Koeppl, 2019) approaches. More recently, Seifner and Sánchez (2023) used neural variational inference (Kingma and Welling, 2013) and neural ODEs (Chen et al., 2018) to infer an implicit distribution over the MJP transition rates. All these variational methods encode the MJP inductive bias into their training objective and, in some cases, into their architecture too.Besides the model of Seifner and Sánchez (2023) , which automatically infers the coarse-grained representation X(t) from D-dimensional, countinuous signals, all the solutions above tackle the MJP inference problem directly in coarse-grained space. Yet below, we also investigate the conformational dynamics of physical systems for which the recorded data lies in a continuous space. To approach such type of problems, we will first need to define a coarse-grained representation of the state space of interest. Fortunately for us, there is a large body of works, within the molecular simulation community, precisely dealing with different methods to obtain such representations, and we refer the reader to e.g. Noé et al. (2020) for a review. McGibbon and Pande (2015) , for example, leveraged one such method to infer the MJP transition rates describing a molecular dynamics simulation via maximum likelihood. Alternatively, researchers have also treated the conformational states in these systems as core sets, and inferred phenomenological MJP rates from them (Schütte et al., 2011) , or modelled the fast intra-state events as diffusion processes, indexed by a hidden MJP, and inferred the latter either via MCMC (Kilic et al., 2021; Köhs et al., 2022) or variational (Horenko et al., 2006; Köhs et al., 2021) methods.In this work we tackle the classical MJP inference problem on coarse-grained space and present, to the best of our knowledge, its first zero-shot solution. 
foundation inference models: 
 In this section we introduce a novel methodology for zero-shot inference of Markov jump processes which frames the inference task as a supervised learning problem. Our main assumption is that the space of realizable MJPs3 , which take values on bounded state spaces that are not too large, is simple enough to be covered by a heuristically constructed synthetic distribution over noisy and discretely observed MJPs. If this assumption were to hold, a model trained to infer the hidden MJPs within a synthetic dataset sampled from this distribution would automatically perform zero-shot inference on any unseen sequence of empirical observations. We do not intend to formally prove this assumption. Rather, we will empirically demonstrate that a model trained in such a way can indeed perform zero-shot inference of MJPs in a variety of cases.Our methodology has two components. First, a data generation model that encodes our believes about the class of realizable MJPs we aim to model. Second, a neural recognition model that maps subsets of the simulated MJP observations onto the initial condition and rate matrix of their target MJPs. We will explore the details of these two components in the following sections.In this subsection we define a broad distribution over possible MJPs, observation times and noise mechanisms, with which we simulate an ensemble of noisy, discretely observed MJPs. Before we start, let us remark that we will slightly abuse notation and denote both probability distributions and their densities with the same symbols. Similarly, we will also denote both random variables and their values with the same symbols.Let us denote the size of the largest state space we include in our ensemble with C, and arrange all transition rates, for every MJPs within the ensemble, into C × C rate matrices. Let us label these matrices with F. We define the probability of recording the noisy sequence x ′ 1 , . . . , x ′ l ∈ X , at the observation times 0 < τ 1 < • • • < τ l < T , with T the observation time horizon, as followsEQUATIONNext, we specify the different components of Eq. 2, starting from the right.Distribution over initial conditions. The distribution p(π 0 |ρ 0 ), with hyperparameter ρ 0 , is defined over the C-simplex, and encodes our beliefs about the initial state (i.e. the preparation) of the system. It enters the master equation as the class probabilities of the categorical distribution over the states of the system, at the start of the process. That is p MJP (x, t = 0) = Cat(π 0 ). We either choose π 0 to be the class probabilities of the stationary distribution of the process, or sample it from a Dirichlet distribution. Appendix B provides the specifics.Distribution over rate matrices. The distribution p rates (F|A, ρ f ) over the rate matrices encodes our beliefs about the class of MJPs we expect to find in practice. We define it to cover MJPs with state spaces whose sizes range from 2 until C, because we want our FIM to be able to handle processes taking values in all those spaces. The distribution is conditioned on the adjacency matrix A, which encodes only connected state spaces (i.e. irreducible embedded Markov chains only), and a hyperparameter ρ f which encodes the range of rate values within the ensemble. Specifically, we define the transition rates asF ij = a ij f ij, where a ij is the corresponding entry of A and f ij is sampled from a set of Beta distributions, with different hyperparameters ρ f . Note that these choices restrict the values of the transition rates within the ensemble to the interval (0, 1) and hence, they restrict the number of resolvable transitions within the time horizon T of the simulation. We refer the reader to Appendix B, where we specify the prior p(A, ρ f ) = p(A)p(ρ f ) and its consequences, as well as give details about the sampling procedure. We also discuss the main limitations of choosing a Beta prior over the transition rates in Section 6.Distribution over observation grids. The distribution p grid (τ 1 , . . . , τ l |ρ τ ), with hyperparameter ρ τ , gives the probability of observing the MJP at the times τ 1 , . . . , τ l , and thus encodes our uncertainty about the recording process. Given that we do not know a priori whether the data will be recorded regularly or irregularly in time, nor we know its recording frequency, we define this distribution to cover both regular and irregular cases, as well as various recording frequencies. Note that the number of observation points on the grid is variable. Please see Appendix B for details.Distribution over noise process. Just as the (instantaneous) solution of the master equation p MJP (x|t, F, π 0 ), the noise distribution p noise (x ′ |x, ρ x ), with hyperparameter ρ x , is defined over the set of metastable states X . Recall that FIM solves the MJP inference problem directly in coarse-grained space. The noise distributions then encodes both, possible measurement errors that propagate through the coarse-grained representation, or noise in the coarse-grained representation itself. We provide details of its implementation in Appendix B.We use the generative model, Eq. 2 above, to generate N MJPs, taking values on state spaces with sizes ranging from 2 to C. We then sample K paths per MJP, with probability p(K), on the interval [0, T ]. The jth instance of the dataset thus consists of K paths and is given byEQUATIONandx ′ jki ∼ p noise (x ′ |X jk (τ jki )) (K,l) (k,i)=(1,1), with τ jk1 , . . . , τ jklK k=1 ∼ p grid (τ 1 , . . . , τ l |ρ τ ),where Gillespie denotes the Gillespie algorithm we use to sample the MJP paths (see Algorithm 1). Note that we make the number of paths (K above) per MJP random, because we do not know a priori how many realizations (i.e. experiments), from the empirical process of interest, will be available at the inference time. We refer the reader to Appendix B for additional details.Figure 2 illustrates the complete data generation process.In this subsection we introduce a neural recognition model that processes a set of K time series of the form {(x ′ k1 , τ k1 ), . . . , (x ′ kl , τ kl )} K k=1 , as generated by the procedure in Eq. 3 above, and estimates the intensity rate matrix F and initial distribution π 0 of the hidden MJP. Practically speaking, we would like the model to be able to infer MJPs from time series with observation times on any scale.To ensure this, we first normalize all observation times to lie on the unit interval, by dividing them by the maximum observation time τ max = max{τ k1 , . . . , τ kl } K k=1 , and then rescale the output of the model accordingly (see Appendix C for details).Let us use ϕ, ψ and Ω to denote feed-forward, sequence processing networks, and attention networks, respectively. Thus ψ can denote e.g. LSTM or Transformer networks, while Ω can denote e.g. a self-attention mechanism. Let us also denote the networks' parameters with θ.We first process each time series with a network ψ 1 to get a set of K embeddings, which we then summarize into a global representation h θ through the attention network Ω 1 . In equations, we writeEQUATIONNext we use the global representation to get an estimate of the intensity rate matrix, which we artificially model as a Gaussian variable with positive mean, and the initial distribution of the hidden MJP as followsEQUATIONwhere the exponential function ensures the positivity of our estimates, and the variance is used to represent the model's uncertainty in the estimation of the rates (Seifner et al., 2024) . The right panel of Figure 2 summarizes the recognition model, and Appendix C provides additional information about the inputs to, outputs of and rescalings done by the model.Training objective. We train the model to maximize the likelihood of its predictions, taking care of the exact zeros (i.e. the missing links) in the data. To witEQUATIONwhere the second term is nothing but the mean-squared error of the predicted rates fij (and its standard deviation) when the corresponding link is missing, and can be understood as a regularizer with weight λ. The latter is a hyperparameter.During training, FIM processes a variable number K of time series, which lies on the interval [K min , K max ]. Similarly, each one of these time series has a variable number l of observation points, which lies on the interval [l min , l max ]. We shall say that FIM needs a bare minimum of K min l min input data points to function. Perhaps unsurprisingly, we have empirically seen that FIM perform bests when processing K max l max data points. Going significantly beyond this number seems nevertheless to decrease the performance of FIM. We invite the reader to check Appendix D for details.Let us define then, for the sake of convenience, the FIM context number c(K, l) = Kl as the number of input points4 FIM makes use of to estimate F and π 0 .The (a priori known) dimension of the process as an integer between 2 and C. If this dimension is unknown, the model returns a C × C rate matrix whose rank might (approximately) be smaller than C, which indicates a hidden state-space of size smaller than C.We recommend users to use the model only within its training range. That is, with up to a maximum of K = 300 paths, and grids up to a maximum of L = 100 points. 
experiments: 
 In this section we test our methodology on five datasets of varying complexity, and corrupted by noise signals of very different nature, whose hidden MJPs are known to take values in state spaces of different sizes. In what follows we use one and the same (pretrained) FIM to infer hidden MJPs from all these datasets, without any parameter fine-tuning. Our FIM was (pre)trained on a dataset of 45K MJPs, defined over state spaces whose sizes range from 2 to 6. A maximum of (K =)300 realizations (paths) per MJP were observed during training, everyone of which spanned a time-horizon T = 10, recorded at a maximum of 100 time points, 1% of which were mislabeled. Given these specifications, FIM is expected to perform best for the context number c(300, 100) during evaluation. Additional GROUND TRUTH 1.00 1.00 1.00 NEURALMJP 1.06 1.17 1.14 FIM 1.11(7) 0.99(8) 0.98(5) Baselines: Depending on the dataset, we compare our findings against the NeuralMJP model of Seifner and Sánchez (2023) , the switching diffusion model (SDiff) of Köhs et al. (2021) , and the discrete-time Markov model (VampNets) of Mardt et al. (2017) .All these baselines are trained on the target datasets.In statistical physics, the ratchet effect refers to the rectification of thermal fluctuations into directed motion to produce work, and goes all the way back to Feynman (Feynman et al., 1965) . Here we consider a simple example thereof, in which a Brownian particle, immersed in a thermal bath at unit temperature, moves on a one-dimensional lattice. The particle is subject to a linear, periodic and asymmetric potential of maximum height 2V that is switched on and off at a constant rate r. The potential has three possible values when is switched on, which correspond to three of the states of the system. The particle jumps among them with rate f on ij . When the potential is switched off, the particle jumps freely with rate f off ij . We can therefore think of the system as a six-state system, as illustrated in Figure 3 . Similar to Roldán and Parrondo (2010) , we now define the transition rates asEQUATIONGiven these specifics, we consider the parameter set (V, r, B) = (1, 1, 1) together with the dataset simulated by Seifner and Sánchez (2023) , which consists of 5000 paths (in coarse-grained space) recorded on an irregular grid of 50 time points. The task is to infer (V, r, B) from these time series. NeuralMJP infers a global distribution over the rate matrices and hence relies on their entire train set, which amounts to about 4500 time series. We therefore report FIM evaluations with context number c(300, 50) on that same train set, averaged over 15 (non-overlapping) batches in Table 1 .The results show that FIM performs on par with (or even better than) NeuralMJP, despite not having been trained on the data. Note in particular that our results are sharply peaked around their mean, indicating that a context of c(300, 50) points only contains enough information to describe the data well. What is more, Table 16 in the Appendix demonstrates that FIM can infer vanishing transition rates as well (see Eq. 6). Now, being able to infer the rate matrix in zero-shot mode allows us to immediately estimate a number of observables of interest without any training. Stationary distributions, relaxation times and mean first-passage times (see Appendix A for their definition), as well as time-dependent moments, can all be computed zero-shot via FIM. For example, we report on the left block of Figure 4 the time-dependent class probabilities (i.e. the master eq. solutions) computed with the FIM-inferred rate matrix (black), against the ground-truth solution (blue). The agreement is very good.Zero-shot estimation of entropy production. The DFR model is interesting because the random switching combined with the asymmetry in the potential make it more likely for the particle to jump towards the right (see Figure 4 ). Indeed, that is the ratchet effect. As a consequence, the system features a stationary distribution with a net current -the so-called non-equilibrium steady state ( Ajdari and Prost, 1992) , which is characterized by a non-vanishing (stochastic) entropy production. The development of (neural) estimators of entropy production is a very active topic of current research (see e.g. Kim et al. (2020) and Otsubo et al. (2022) ). Given that the entropy production can be written down in closed form as a function of both the rate matrix and the master eq. solution (see e.g. Seifert ( 2012)), we can readily use FIM to estimate it.Figure 4 displays the total entropy production computed with FIM for a set of different potentials. The results are averaged over 15 FIM evaluations with c(300, 50) and are again in very good agreement with the ground truth. It is noteworthy that FIM, trained on our heuristically constructed dataset, captures well a continuous set of MJPs. That is, we evaluate one and the same FIM over different datasets, each sampled from a DFR model with a different potential value. In sharp contrast, state-ofthe-art models need to be retrained for every new potential value (Kim et al., 2020) .Zero-shot simulation of the DFR process. Inferring the rate matrix and initial condition of a MJP process entails that one can also sample from it. Our FIM can thus be used as a zero-shot generative model for MJPs. However, to test the quality of said MJP realizations wrt. some target MJP, we need a distance between the two. Here we propose to use the Hellinger distance (Le Cam and Yang, 2000) to first estimate the divergence between a sequence of (local) histogram pairs, recorded at a given set of observation times, and then average the local estimates along time. Appendix F.1 empirically demonstrates that this pragmatically defined MJP distance is sensible.Table 2 reports the time-averaged Hellinger distance between 1000 (ground-truth) DFR paths and 1000 paths sampled from (the MJPs inferred by) NeuralMJP and FIM. We repeat this calculation 100 times, for 1000 newly sampled paths from NeuralMJP and FIM, but the same 1000 target paths, to compute the mean values and error bars in the Table . The results show that the zero-shot DFR simulation obtained through FIM is on par with the NeuralMJP-based simulation, wrt. the ground truth.In this section we study the conformational dynamics of the viral ion channel Kcv MT325 , which exhibits three metastable states (Gazzarrini et al., 2006) . Specifically, we analyse the ion flow across the membrane as the system jumps between its metastable configurations. This ion flow was recorded at a frequency of 5kHz over one second. Figure 1 shows one snapshot of these recordings, which were made available to us via private communication (see the Acknowledgements). Our goal is to infer physical observables -like the stationary distribution and mean first-passage times -of the conformational dynamics, and to compare our findings against the SDiff model of Köhs et al. (2021) and NeuralMJP.The recordings live in real space, which means that we first need to obtain a coarse-grained representation (CGR) from them, before we can apply FIM. Here we consider two CGRs: the CGR inferred by NeuralMJP and a naive CGR obtained with a Gaussian Mixture Model (GMM). Given that we only have 5000 observations available, we make use of a single FIM evaluation with context number c(50, 100). We infer two FIM rate matrices, one per each CGR, which we label as FIM-NMJP and FIM-GMM.Table 3 contains the inferred stationary distributions from all models and evidences that a single FIM evaluation is enough to unveil the long-time asymptotics of the process. Similarly, Table 15 in the Appendix, which contains the inferred mean-first passage times, demonstrates that FIM makes the same inference about the short-term dynamics of the process as do SDiff and NeuralMJP. See Appendix F for additional results.Zero-shot simulation of switching ion channel process. Just as we did with the DFR process, we can use FIM to simulate the switching ion channel process in coarse-grained space. Since only paths on the same CG space can be compared, we evaluate NeuralMJP against FIM-NMJP. To construct the target distribution, we leverage another 30 seconds of measurements, which amount to 150K observations that have not been seen by any of the models. The results in Table 2 indicate that our zero-shot simulations is statistically closer to the ground-truth process than the NeuralMJP simulation.Alanine dipeptide is 22-atom molecule widely used as benchmark in molecular dynamics simulation studies. Its popularity stems from the fact that the heavy-atom dynamics, which jumps between six metastable states, can be fully described in terms of the dihedral (torsional) angles ψ and ϕ (see e.g. Mironov et al. (2019) for details).We examine an all-atom ADP simulation of 1 microsecond, which was made available to us via private communication (see the Acknowledgements below), and compare against both, the VampNets model of Mardt et al. (2017) and NeuralMJP. The data consists of the values taken by the dihedral angles as time evolves and thus needs to be mapped onto some coarse-grained space. We again make use of NeuralMJP to obtain a CGR. We then use FIM with context number c(300, 100) to process 32 100-point time windows of the simulation and compute an average rate matrix. Note that this is the optimal context number of our pretrained model. Table 4 (and Appendix F.2) confirms that, once again, FIM can infer the same physical properties from the ADP simulation as the baselines.Zero-shot simulation of the alanine dipeptide. Simulations in coarse-grained space for molecular dynamics is a high-interest research direction (Husic et al., 2020) . Here we demonstrate that FIM can be used to simulate the ADP process in zero-shot mode. Indeed, Table 2 reports the distance from both NeuralMJP and FIM to a target ADP process, computed from 200 paths with 100 observations each. Once more, FIM performs comparable to NeuralMJP.Finally, we consider two additional systems that feature jumps between two metastable states: a simple protein folding model and a two-mode switching system. We invite the reader to check out Appendix F.5 and F.6 for the details. That being said,  
conclusions: 
 In this work we introduced a novel methodology for zero-shot inference of Markov jump processes and its Foundation Inference Model (FIM). We empirically demonstrated that one and the same FIM can be used to estimate stationary distributions, relaxation times, mean first-passage times, time-dependent moments and thermodynamic quantities (i.e. the entropy production) from noisy and discretely observed MJPs, taking values in state spaces of different dimensionalities, all in zero-shot mode. To the best of our knowledge, FIM is also the first zero-shot generative model for MJPs.Future work shall involve extending our methodology to Birth and Death processes, as well as considering more complex (prior) transition rate distributions. See our discussion on Limitations in the next section, for details. 
limitations: 
 The main limitations of our methodology clearly involve our synthetic distribution. Evaluating FIM on empirical datasets whose distribution significantly deviates from our synthetic distribution will, inevitably, yield poor estimates. Consider Figure 4 (right), for example. The performance of FIM quickly deteriorates for V ≥ 3, for which the ratio between the largest and smallest rates gets larger than about three orders of magnitude. These cases are unlikely under our prior Beta distributions, and hence effectively lie outside of our synthetic distribution.More generally, the MJP dynamics underlying phenomena that feature long-lived, metastable states, ultimately depends on the shape of the energy landscape characterizing the set X , inasmuch as the transition rates between metastable states i and j (f ij in our notation) are characterized by the depth of the energy traps (that is, the height of the barrier between them).In equations, we writeEQUATIONwhere E j is the jth trap depth, and T is the temperature of the system. Therefore, the distribution over energy traps determines the distribution over transition rates.  
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: Our paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: We have credited the owners of the evaluation data and referenced the related work on which this project has been built on.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 13. New Assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [Yes] Justification: The new asset of this paper are the code and the models which are well documented.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: In this section we provide some brief background on MJPs and describe how physical quantities such as the stationary distributions, relaxation times and mean first passage times can be computed from the intensity matrix. Additionally, we mention how trajectories for MJPs can be sampled using the Gillespie algorithm.Markov jump processes are stochastic models used to describe systems that transition between states at random times. These processes are characterized by the Markov property where the future state depends only on the current state, not on the sequence of events that preceded it.A continuous-time MJP X(t) has right-continuous, piecewise-constant paths and takes values in a countable state space X over a time interval [0, T ]. The instantaneous probability rate of transitioning from state x ′ to x is defined asEQUATIONwhere p MJP (x, t|x ′ , t ′ ) denotes the transition probability.The evolution of the state probabilities p MJP (x, t) is governed by the master equationEQUATION)For homogeneous MJPs with time-independent transition rates, the master equation in matrix form isEQUATIONwith the solution given by the matrix exponentialEQUATIONThe stationary distribution p * MJP of a homogeneous MJP is a probability distribution over the state space X that satisfies the condition p * MJP • F = 0. This implies that the stationary distribution is a left eigenvector of the rate matrix corresponding to the eigenvalue 0.The relaxation time of a homogeneous MJP is determined by its non-zero eigenvalues λ 2 , λ 3 , . . . , λ |X | . These eigenvalues define the time scales of the process:|Re(λ 2 )| -1 , |Re(λ 3 )| -1 , . . . , |Re(λ |X | )| -1 .These time scales are indicative of the exponential rates of decay toward the stationary distribution. The relaxation time, which is the longest of these time scales, dominates the long-term convergence behavior. If the eigenvalue corresponding to the relaxation time has a non-zero imaginary part, then this means that the system does not converge into a fixed stationary distribution but that it instead ends in a periodic oscillation.For an MJP starting in a state i ∈ X , the first-passage time to another state j ∈ X is defined as the earliest time t at which the MJP reaches state j, given it started in state i. The mean first-passage time (MFPT) τ ij is the expected value of this time. For a finite state, time-homogeneous MJP, the MFPTs can be determined by solving a series of linear equations for each state j, distinct from i, with the initial condition that τ ii = 0EQUATIONA.5 The Gillespie Algorithm for Continuous-Time Markov Jump ProcessesThe Gillespie algorithm (Gillespie, 1977) is a stochastic simulation algorithm used to generate trajectories of Markov jump processes in continuous time. The algorithm proceeds as follows:Algorithm 1 Gillespie Algorithm for Markov Jump Processes 1: INPUT: The intensity matrix F, the initial state distribution π 0 , the starting time t 0 and the end time t end 2: Initialize the time t to the starting time t 0 3: Initialize the system's state s to an initial state s 0 ∼ π 0 4: While t < t end do 5:Calculate the intensity λ = -1/F ss from state s 6:Sample the time τ to the next event from an exponential distribution with rate λ Sample the next state s ′ from the distribution defined by p 12:Update the system's state s ← s ′ 13: Record the state s and time t 14: End while 15: OUTPUT: The trajectory of states and timesThis section is a continuation of section 3.1 and provides more details on the generation of our synthetic training dataset. Additionally, we provide some statistics about the dataset distribution.In this subsection we give additional details about our data generation mechanism.Distribution over rate matrices. Our data generation procedure starts by sampling the entries f ij of the intensity matrix from the following beta distributionsp(f ij |ρ f ) = Beta(ρ f = (α, β)), with p(α) = Uniform({1, 2})and p(β) = Uniform({1, 3, 5, 10}). ( 14)Both these discrete uniform distribution define the prior p(ρ f ) = p(α)p(β).The choices for α and β were made heuristically, to obtain reasonable (i.e. varied) distributions over the number of jumps (see e.g. Figure 5 ). We remark that we fixed this set of training distributions before evaluating the model on the evaluation sets, in order to prevent us from introducing unwanted biases into the distribution hyperparameters by optimizing on the evaluation set.Next we define the prior over the adjacency matrix asEQUATIONwhere δ(•) labels the Dirac delta distribution and J denotes the matrix for which all off-diagonal entries are 1 and the diagonal ones are 0. Furthermore p Erdös-Rényi labels the Erdös-Rényi model (Erdös and Rényi, 1959) , for which each link is defined via an independent Bernoulli variable, with some fixed, global probability p, here set to 1 2 . Equation 15indicates that (in average) 50 percent of our state networks are fully connected, whether the other 50 percent are not.Our motivation for this prior is that it often happens in real world processes that the intensity matrices are not fully connected. Let us remark, however, that we only accept the Erdös-Rényi sample if the corresponding graph is connected -that is, if the system cannot get stuck into a single state. Both these distributions implicitly define p rates (F|A, ρ f ), for F ij = a ij f ij . Remark on generalization beyond prior rate distribution. We remark that while all entries of the intensity matrix seen during training lie on the interval [0, 1], the model can still predict intensities outside this interval. We empirically demonstrated that this in indeed the case on the widely different target sets of the experimental section, in the main text. The reason behind this is that we normalize the maximum time among all input paths to be 1, and rescale the predicted intensities accordingly.Ultimately, what matters is the difference among the rates (and therefore among the observation times) within the target time series. Our approach for sampling intensity matrices resulted in a vast variety of different processes.The distribution of the number of jumps per trajectory is shown in figure 5 and that of relaxation times is shown in figure 6 .Distribution over initial conditions. We choose half of our initial distributions in our synthetic ensemble to be the stationary distribution of the MJP p * MJP . The motivation for this is that it often happens that real life experiments produce very long observations of a system in equilibrium. The second half of our initial distributions π 0 are randomly sampled from a Dirichlet distribution Dir(ρ 0 ), where we heuristically choose ρ 0 = 50. In equations, we writeEQUATION)Distribution over observation grids. In practice, the exact jump (i.e. transition) times are not known. We therefore first generate observations of the state of the system on a regular grid with a maximum of L = 100 points. We then randomly mask out some observations from this fixed regular grid, in order to make the model grid independent. Half of our (subsampled) observation grids are chosen to be regular, i.e. they are strided with strides ∈ {1, 2, 3, 4}. The other half are chosen to be irregular, through a Bernoulli filter (or mask) with ρ survival ∈ {1/4, 1/2} applied to the base (L = 100) grid.Distribution over noise process. Because real world data is often noisy we also add noise to the labels. If a state observation is selected to be mislabeled, the new label is randomly chosen from a uniform distribution over all states. We investigate two different configurations in this project, one with 1% label noise (ρ x = 0.01) and one with 10% label noise (ρ x = 0.1).MJP simulation. We sample the jumps between different states with an algorithm due to (Gillespie, 1977) (see A.5). We sample jumps between times 0 and 10 because almost all of our processes are in equilibrium by then (see figure 6 ). Figure 6 : Distributions of the relaxation times. We also report the percentage of processes that converge into an oscillating distribution (OP) and the percentage of processes that have a relaxation time which is larger than the maximum sampling time (NCP) of our training data (given by t end = 10). The figures are based on 1000 processes.Training Dataset Size The synthetic dataset on which our models were trained consists of 25k six-state processes, and 5k processes of 2-5 states, resulting in a total size of 45k processes. For each of these processes we sampled 300 paths.Distribution over the number of MJP paths p(K). While we generate the data with 300 paths per process, we want to ensure that the model is able to handle datasets with less than 300 paths. For this reason, we shuffle the training data at the beginning of every epoch and distribute it into batches with path counts 1, 11, 21, . . . , 300. We found that such a static selection of the path counts is better than a random selection, because a random selection leads to oscillating loss functions (because the model obviously gets a larger loss for samples with fewer paths), and thus training instabilities. Since we do not always select all paths per process but instead select a random subset of them, the data that the model processes changes during every epoch, which helps in reducing overfitting.In this section, we give details about the inputs to and outputs of our pretrained recognition model. We also comment on the internal rescalings done by the model, in order to be able to infer MJPs from time series with observation times of any scale.The model takes as input three parameters:1. The observation grids (shape: [num_paths K, grid_size L]): The observation times {τ k1 , . . . , τ kl } K k=1 , padded to the maximum length L. 2. The observation values (shape: [num_paths K, grid_size L]): The noisy observation values {x ′ k1 , . . . , x ′ kl } K k=1 padded to the maximum length L. Note that these values are integers lying on the discrete set {0, 1, . . . , C -1}.Internally, the model does the following:1. It computes the maximum observation time:EQUATION2. It normalizes the observation times between 0 and 1:EQUATION3. It computes the inter-event times ∆τ ki = τ k,i+1 -τ ki , for k = 1, . . . , K.4. It transforms the observation values to one-hot-encodings.5. It predicts the (normalized) off-diagonal elements of the intensity matrix and variance matrix as well as the initial distribution (here we are working with the maximum supported dimension, that is C).6. It rescales back the estimates to the original time scale:EQUATIONEQUATIONNote that, as we empirically demonstrated in the paper, this rescaling procedure allows us to work with real-world MJPs of arbitrary time scales. For example, the time scales for the switching ion channel dataset were more than 500 times smaller than the time scales in our training dataset.We now elaborate on how the model can deal with processes whose state spaces have sizes c < C.We arranged all the target rate matrices F within our training dataset, for MJPs with state spaces of size c < C, to be the leftmost block diagonal c × c matrix within a C × C matrix of zeros, so that the redundant matrix elements are always zero. As can be read from equation ( 6) of the main text, we train FIM to predict zeros for those redundant matrix elements.In practice, however, our trained FIM does not exactly predict zeros for those redundant matrix elements. In our experiments, the user knows a priori the number of states c of the hidden process, so we explicitly set the redundant matrix elements to zero, and only then compute the corrected diagonal (i.e. the normalization) of the output rate matrix.We afterwards select the c × c entries of the predicted intensity matrix and variance matrix as well as the first c entries of the predicted initial distribution.We refer the reader to our library for additional details.The output of the model consists of three parameters:1. The intensity matrix F (shape: [c, c]).2. The variance matrix Var F (shape: [c,c]).3. The initial distribution π 0 (shape [c]).In this section we provide more details about the architecture of our models and the hyperparameters.Path encoder ψ 1 . We evaluated two different approaches for the path encoder ψ 1 . The first approach utilizes a bidirectional LSTM (Hochreiter and Schmidhuber, 1997) as ψ 1 , while the second approach employs a transformer (Vaswani et al., 2017) for ψ 1 . The time series embeddings are denoted by h kθ (see Equation 4). The input to the encoder ψ 1 is (x ′ k1 , τ k1 , . . . , x ′ kl , τ kl ), where τ kl = [τ kl , δ kl ], δ kl = τ kl -τ (k-1)l , and x kl ∈ {0, 1} C is the one-hot encoding of the system's state.Path attention network Ω 1 . We tested two approaches. The first approach uses classical selfattention Vaswani et al. (2017) and selects the last embedding. For the second approach we used an approach we denote as learnable query attention which is equivalent to classical multi-head attention with the exception that we do not compute the query based on the input, but instead make it a learnable parameter, i.e.,EQUATION)EQUATIONwhereH 1:K ∈ R K×d model denotes a concatenation of h 1 , . . . , h K , W K i , W V i ∈ R d model ×d k and Q i ∈ R q×d kis the learnable query matrix. The output dimension of the learnable query attention is therefore independent of the number of input tokens.Hyperparameter tuning: Hyperparameters were tuned using a grid search method. The optimizer utilized was AdamW (Loshchilov and Hutter, 2017) , with a learning rate and weight decay both set at 1e -4 . A batch size of 128 was used. During the grid search, we experimented with the hidden size of the path encoder ([64, 128, 256, 512] ), the hidden size of the path attention network ([128, 256] ), and various MLP architectures for ϕ 1 , ϕ 2 , and ϕ 3 ([[32, 32] , [128, 128] Training procedure: All models were trained on two A100 80Gb GPUs for approximately 500 epochs or approximately 2.5 days on average per model. Early stopping was employed as the stopping criterion. The models were trained by maximizing the likelihood.The final models (FIM-MJP 1% Noise and FIM-MJP 10% Noise) have the following hyperparameters: Path encoder -hidden_size(ψ 1 ) = 256 (the final models used a BiLSTM); Path attention network -Ω 1 : q = 16, d k = 128 (the final models used the learnable query approach); ϕ 1 , ϕ 2 , ϕ 3 = [128, 128] .Pretrained models: Our pretrained models are also available online5 .In this section, we study the performance of the models with different architectures. Additionally, we study the behavior of the performance of the models with respect to varying numbers of states and varying number of paths.If the evaluation set is larger than the optimal context number c(K max , l max ), we split the evaluation set into batches and give these to the model independently (because the model does not work well to give the model more paths than during training, see table 8 ). Afterwards, we compute the mean of the predictions among the batches and report the mean RMSE of the intensity entries (if the ground-truth is available). This makes it easier to compare our model against previous works which have also used the full dataset to make predictions. Interestingly, we find that the RMSE of this averaged prediction is often significantly better than the mean RMSE among the batches. For example for the DFR dataset the RMSE of the averaged prediction is 0.0617, while the average RMSE of the batches is 0.122. If the dataset has been split into multiple batches, we report the RMSE together with the standard deviation of the RMSE among the batches. The reported confidence is the mean predicted variance of the model (recall that we are using Gaussian log-likelihood during training).The ablation study presented in Table 5 evaluates the impact of different model features on the performance by comparing various combinations of architectures and attention mechanisms with varying numbers of paths, and their corresponding RMSE values. The study examines models using a BiLSTM or Transformer, with and without self-attention and learnable query attention, across 1, 100, and 300 paths. The results indicate that increasing the number of paths consistently reduces RMSE (see section E.4 for more details), demonstrating the benefit of considering more paths during training. Specifically, using a BiLSTM with learnable query attention achieves an RMSE of 0.193 ± 0.031 with a single path, significantly improving to 0.048 ± 0.011 with 100 paths, and further to 0.0457 ± 0.0 with 300 paths. Similarly, a Transformer with learnable query attention shows an RMSE of 0.196 ± 0.031 for a single path, 0.049 ± 0.011 for 100 paths, and 0.0458 ± 0.0 for 300 paths. The inclusion of self-attention in the Transformer models slightly improves performance, with the best RMSE of 0.0459 ± 0.0 achieved when both self-attention and learnable query attention are used with 300 paths. In this case since many of the processes contain one path it is beneficial to use the learnable query attention over the standard self-attention mechanism. (64, 128, 256, and 512) . The RMSE remains relatively stable across different hidden sizes of Ω 1 , with slight variations observed depending on the hidden size of ψ 1 . Overall, the plots highlight that some components, such as ψ 1 and ϕ 1 , are more sensitive to changes in hyperparameters, emphasizing the importance of selecting appropriate hyperparameters to optimize model performance. ). The third plot shows minimal RMSE impact from ϕ 2 architecture size. The fourth plot shows RMSE stability across different Ω 1 hidden sizes, with slight variations based on ψ 1 . This highlights the importance of tuning ψ 1 and ϕ 1 for optimal performance. 1% Noise Data 10% Noise Data FIM-MJP 1% 0.046 0.199 FIM-MJP 10% 0.096 0.087Table 6 : Performance of FIM-MJP 1% and FIM-MJP 10% on synthetic datasets with different noise levels. We use a weighted average among the datasets with different numbers of states to compute a final RMSE.Table 6 compares the performance of two models, FIM-MJP 1% and FIM-MJP 10%, on synthetic datasets with noise levels of 1% and 10%, measured in terms of RMSE. For datasets with 1% noise, the FIM-MJP 1% model achieves an RMSE of 0.046, indicating good performance, but its RMSE increases significantly to 0.199 on 10% noise data, showing decreased performance with higher noise. Conversely, the FIM-MJP 10% model, trained with 10% noise data, has an RMSE of 0.096 on 1% noise data, higher than the FIM-MJP 1% model on the same data, but achieves a lower RMSE of 0.087 on 10% noise data, demonstrating better performance under high noise conditions. This indicates that the FIM-MJP 10% model is more robust to noise, maintaining consistent performance across varying noise levels, while the FIM-MJP 1% model excels in low noise environments but struggles with higher noise. The results highlight the importance of training with appropriate noise levels to ensure robust model performance across different noise conditions.We compare the performance of our models on processes with varying number of states. Note that our model always outputs a 6 × 6 dimensional intensity matrix. However, in these experiments we only use the rows and columns that correspond to the lower-dimensional process. This improves the comparability between different dimensions as lower-dimensional processes obviously have many zero-entries in their intensity matrix which would make it easier for the model to achieve a good RMSE score.It can be seen in Table 7 that the multi-state-model performs well among all different dimensions. As expected, lower-dimensional processes seem to be easier for the model. Additionally, Table 7 shows the performance of a model which has only been trained on six-state processes. The performance of this native six-state-model for six number of states is very similar to the multi-state-model which shows that having more states during training does not reduce the single-state performance. As expected, the performance of the six-state model on processes with lower numbers of states is significantly worse, but still better than random.One of the advantages of our model architecture is that it can handle arbitrary number of paths.We therefore use our model that was trained on at maximum 300 paths and assess its performance 8 . When being inside the training range, the performance and the confidence of the model goes down as the model is given fewer paths per evaluation, which is to be expected. Interestingly, the performance of the learnable-query (LQ) model peaks at 500 paths instead of at 300, which was the maximum training range. One possible explanation for this might be that we are still close enough to the training range while being able to use the full data (note that the dataset contains 5000 paths which is not divisible by 300, so we have to leave some of the data out). Going too far beyond the training range does however not work well, for example processing all 5000 paths at once leads to very poor performance, although the model (falsely) become very confident. Another insight from this experiment is that the self-attention (SA) architecture behaves significantly worse when going beyond the maximum number of paths that was seen during training. This is another reason why we chose the (LQ) architecture over the (SA) architecture for the final version of our model. This section contains more of our results which did not fit into the main text. We begin this section by providing more details on the Hellinger distance which we used as a metric to assess the performance of our models. Afterwards, we provide more results and background on the ADP, ion channel and DFR datasets. Additionally, we introduce two two-state MJPs, given by the protein folding datasets (F.5) and the two-mode switching system (F.6), which we use to evaluate our models and to compare it against previous works.Real-world empirical datasets of MJPs provide no knowledge of a ground truth solution. For this reason we present a new metric that can be used to compare the performance of the inference of various models based on only the empirical data. Our metric of choice is the Hellinger distance which is a measure of the dissimilarity between two probability distributions. Given two discrete probability distributions P = (p 1 , . . . , p k ) and Q = (q 1 , . . . , q k ), the Hellinger distance is defined asEQUATION)For our empirical cases, the class probabilities of the discrete probability distributions are not known explicitly. We therefore approximate them by using the empirical distributions, given by the (normalized) histograms of the observed states at the observation grids.We test this approach on the DFR process by first sampling a specified number of paths for the potential V = 1 using the Gillespie algorithm, which we then consider as the target distribution.Counting states among the different paths then yields histograms of the states for every time step. We repeat the same procedure for different choices of V . Afterwards we compute the Hellinger distance between the newly sampled histogram and the target distribution for every time step. Figure 8 shows that the distance indeed goes down as we approach the target distribution, which provides heuristic evidence of the effectiveness of our metric. The Hellinger distances for various models are shown in Table 2 and Table 9 .As one can see, FIM-MJP performs as well (and sometimes better) as the current state-of-the-art model NeuralMJP. The plot shows the Hellinger distance to a target dataset that was sampled from a DFR with V = 1 on a grid of 50 points between 0 and 2.5. The means and standard deviations were computed by sampling 100 histograms per dataset. As expected, the distance decreases as the voltage gets closer to the voltage of the target dataset. We also remark that the scale of the distances gets smaller as one takes more paths into account and converge to the distance of the solutions of the master equation. 9 : Comparison of the time-average Hellinger distances for various models. We used the same labels as NeuralMJP to make the results comparable. The errors are the standard deviation among 100 sampled histograms. The target datasets contain 200 paths for ADP, 1500 paths for Ion Channel, 2000 paths for Protein Folding and 1000 paths for the DFR. The distances are reported in a scale 1e-2. We remark that the high variance of the distances on the Protein Folding dataset is caused by the models performing basically perfect predictions, which causes the oscillations to be noise. We verified this claim by confirming that the distances of the predictions of the models are as small as the distance of the target dataset to additional simulated data.We use the dataset of Husic et al. (2020) , which models the conformal dynamics of ADP, for evaluating our model. This dataset was provided to us via private communication. The dataset consists of 9800 paths on grids of size 100 and has the sines and cosines of the Ramachandran angles as features:sin ψ, cos ψ, sin ϕ and cos ϕ. We use KMeans to classify the data into states. The reason why we did not choose GMM as for the other datasets is that we could initialize KMeans with hand-selected values to try to achieve a similar classification like those learned by NeuralMJP (Seifner and Sánchez, 2023) , see Figure 9 . Still, the classification is very different and thus also leads to very different results (see Table 10 ). We use 9600 paths to evaluate our models. Our results are shown in Table 10 . Table 11 reports the stationary distributions and compares them to previous works, while Table 12 reports the ordered time scales. Table 15: Mean first-passage times of the predictions of various models on the Switching Ion Channel dataset. We compare against (Köhs et al., 2021) and NeuralMJP (Seifner and Sánchez, 2023) . Entry j in row i is mean first-passage time of transition i → j of the corresponding model.We use the same datasets that were used in Seifner and Sánchez (2023) , which contains 5000 paths on grids of size 50 that lie between times 0 and 2.5. This dataset was provided to us via private communication. We used 4500 paths to evaluate our model. The predicted intensity matrices for the DFR and the ground truth are shown in table 16 . 0 -3 1 1 0 1 0 1 -3 1 0 0 1 1 1 -3        FIM-MJP 1% Noise        -1.dx(t) = -∇U (x(t)) + √ 2dW (t) ,with the potential U (x) being dependent solely on the norm r(x) = |x| as follows:U (x) = -2.5[r(x) -3] 2 , if r(x) < 3 0.5[r(x) -3] 3 -[r(x) -3] 2 , if r(x) ≥ 3This model exhibits bistability in the norm r(x), encapsulating two states akin to the folded and unfolded conformations of a protein.We use the dataset of Seifner and Sánchez (2023) and apply a Gaussian-Mixture-Model to classify the dataset into two states. The decision boundary of the classifier seems to be based on the absolute absolute value of the radius, namely the classifier seems to classify all states with a radius smaller than approximately 2 into the lower state (see figure 11 ).Seifner and Sánchez (2023) generated 1000 trajectories, each with 100 steps after a 1000-step burn-in period. We used 900 paths to evaluate our model. The results are shown in table 17. Table 18 compares the stationary distributions obtained from our models to the ones from Mardt et al. (2017) and Seifner and Sánchez (2023) . with parameters α 1 = α 2 = 1.5, β 1 = -1, and β 2 = 1. For a concise overview of the generation process, the reader is directed to (Köhs et al., 2021) for comprehensive details. We use the same dataset that was generated in (Seifner and Sánchez, 2023) using the code of (Köhs et al., 2021) • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).https://github.com/cvejoski/OpenFIM 38th Conference on Neural Information Processing Systems (NeurIPS2024).We name our model foundation model because it aligns with the definition proposed byBommasani et al. (2021). Indeed, they define foundation models as any model that is trained on broad data (generally using self-supervision at scale) that can be adapted to a wide range of downstream tasks.By realizable MJPs we mean here MJPs that can be inferred from physical processes, given the typical experimental constraints, like e.g. temporal or spatial resolution.We can think about it as the context length in large language models.https://github.com/cvejoski/OpenFIMhttps://github.com/cvejoski/OpenFIM 


Decision: The reviewers agree that this is an interesting paper that deserves to be published. A limitation is that this work uses empirical rather than theoretical investigations to demonstrate the usefulness of the proposed idea, and the paper would be strengthened if this experiments more clearly demonstrated the limitations of this methodology. The authors might want to think about this (as well as other reviewer comments) when they prepare the final version of this paper.

Reviews: review: paper_summary: paper_summary: This work presents a framework for amortizing inference on Markov Jump Processes by learning a foundational model in a supervised fashion from synthetic data. Once learned, the "foundation model" is shown to be successful at zero-shot inference in MJP across a range of domains, out-performing SOTA models that are fine-tuned to the target datasets.
strength_weakness: strength_weakness: Discussion of prior literature could be improved. Some relevant work is on amortized inference in other types of Bayesian models. E.g. "Inference Networks for Sequential Monte Carlo in Graphical Models" learns a supervised model on synthetic data for zero shot inference on Bayesian regression coefficients.
questions: questions: You mention a main limitation is that generalization is poor outside of the family of models assumed in the synthetic training data. What happens to generalization as you scale up the NN model and widen the family of MJPs in the training data (like for LLMs)?
limitations: limitations: Limitations are adequately discussed.
paper_summary: paper_summary: This study describes a foundation model for a specific stochastic process, called Markov Jump Process (MJP). The foundation model, called FIM, is trained with a large number of synthetically generated MJP data set. It is shown that the pre-trained FIM can make a zero-shot inference. The capability of FIM is demonstrated by using a few data set.
strength_weakness: strength_weakness: The numerical experiments are not strong enough. Numerical experiments should achieve either 1) comparing the performance of the proposed model against a cohort of the state-of-the-art methods, or 2) to systematically demonstrate the capability and limitations of the proposed model and investigate the behavior of the model under a range of important conditions. However, the manuscript does not achieve both. Based on the Appendix, I believe that the authors would have performed an extensive number of studies. A more detailed analysis of the model behaviors with respect to the range of parameters will be appreciated.
questions: questions: 1. It is not clearly explained how the data are scaled. Since the training data is generated from equations, it should be fine. But for an inference the input data should be correctly scaled. How is the inference data scaled? Are there any effects of the scaling on the inference result?

2. Similar to question 1, how the time scale is determined. The model predicts "rate", which depends on the time scale of the data. What method do you use for the inference time? Is the time scaling method general enough to be applied for any real-life MJP?

3. A more detailed explanation about the number of states is required. The model is trained with the maximum number of states of six. What happens when the number of states is less than 6? What do $A$ and $F$ look like? Do the redundant elements of the matrices become exactly zero? If not, how do you normalize the transition probability?

4. Similar to 3, it is interesting to see how the model behaves if the number of states is larger than the maximum number of states in the training, as in real-life problems it is difficult to know or restrict the number of states.

5. During the training, when the number of states is less than six, how do you compute $Var\hat{f}_{ij}$? Note that $Var\hat{f}_{ij}$ is used in the denominator and in a log. 

6. How accurate is the prediction of the initial state, $\pi_0$? 

7. The prediction requires a prior distribution for $A$. What are the effects of $A$ on the prediction? What happens if the data is generated from a MJP with a different distribution of $A$?

8. The inference of FIM requires a fairly large amount of input data. A conventional approach is to fit a state-space model with the same amount of data to learn the MJP parameters and make an inference. How does FIM compare against the standard ML approach?
limitations: limitations: NA
paper_summary: paper_summary: The authors propose a foundation model for Markov jump processes (MJPs) that is trained on sequences drawn from synthetic MJPs to predict the corresponding rate matrix and initial state distribution. High-level arguments are presented for the justification of model trained in such a way to be able to generalize to unseen sequences drawn from different processes. Extensive experimental results are shown demonstrating the model's performance after initial training on synthetic data and without any further fine-tuning.
strength_weakness: strength_weakness: As pointed out in the paper, the model naturally has weaknesses generalizing to data that lies outside of the initially exposed distribution that it was exposed to in training, such as on sequences with very high rates or with states that lie outside the initial support.

Outside of this, I could not find any other weaknesses of the paper; I found it to be very good in general.
questions: questions: I personally believe the success demonstrated by the foundation model with how well it generalizes being due to the simplicity (not a bad thing!) of the data, i.e., single dimensional across time and is Markov. Would you agree with this or am I missing something?

Additionally, I was wondering if "zero-shot" is really appropriate to use here. Typically, zero-shot is used when trying to produce predictions when there is no labeled data to learn from, i.e., predicting "dogs" vs. "cats" with an image model that has seen neither during training. This is commonly achieved by providing some side-information, such as a description of a class. Here the model does indeed not see any "labeled" rate matrices; however, it does see many different sequences realized from this matrix. I see it as less zero-shot and moreso that the inference problem has been amortized. By this, I mean that other methods directly learn a MJP whereas this predicts the parameters to one given a batch of sequences. I am not saying any of this to detract from the paper, but rather just to dial down how to actually think about what it is doing and the language used to describe it. Any thoughts on this would be great to hear.
limitations: limitations: The authors adequately describe the limitations.
paper_summary: paper_summary: The paper explores the possibility of using transformers to directly infer the parameters of a Markov jump process (MJPs) from a noisy time data set obtained from a time-dependent process. They first train the machine to correctly predict the model parameters by using synthetic datasets containing data obtained with different realizations of MJPs defined with different parameter sets. Later, they show that the models trained in this way are nevertheless able to derive the parameters of a real simulated or experimental process, even though they were not trained on this data. They also show that one-shot inference is as good as state-of-the-art methods when trained on the target datasets.
strength_weakness: strength_weakness: Despite the originality of the work, there are some aspects that remain unclear to me. I find that important details about how to apply the trained models to real data in practice are missing. I also think the limitations are not properly explored.
questions: questions: Do they make their predictions on the basis of coarse-grained data? If so, how are these representations obtained? If not, how is the data standardized so that they can transfer the models?

I did not understand if the adjacency matrix is learned or not, and I am a bit confused about the information given to the machine in real datasets if not.

I find that the limitations of this approach are not properly explored or discussed in the text. Are the processes used to generate the synthetic datasets realistic in a general setup? Are they able to describe complex processes like the data generated with the trap model? or just Gaussian processes such as a Brownian motion?
limitations: limitations: I discused my doubts about the limitations above.

Comments: title: Response to Rebuttal
comment: I thank the authors for their thoughtful response to my review and am glad to keep my score as accept. I look forward to future work in this area building upon this paper.
We would like to thank the reviewer again for their comment and respond that we will include the discussion above into the Limitations section, as suggested.
We would like to thank the reviewer for their comment and respond that what the reviewer asks is precisely what we did.

Our three baselines (NeuralMJP, SDiff, and VampNets) were trained in the conventional way. That is, our three baseline models were trained on the target datasets. We mention this in our abstract and again on line 82. See also the discussion on lines 234, 237, and 238.

We realize, however, that we may not have been that clear about this. Therefore, to make this point more explicit, we will add a line in the Experiments section to clarify that all baselines are trained on the target datasets.
I want to thank again the authors for the detailed explanation. That's what I expected. I think the issue should be discussed in the limitations section. Trap like dynamics are relevant in the context of disordered proteins, as the random energy model presents trap dynamics.

I will update my score to accept.
I would like to thank the author for the clarifications. My last comment is a question regarding a model trained for the target distribution versus a foundation model trained with a large number of data. In the conventional way, you train a MJP  model for a particular data set and make an inference for the data set with the same underlying distribution. In FIM, you train FIM using a large collection of data set and make an inference on a new data set. It would be interesting to compare the accuracy of these two different approaches. It will be surprising if FIM outperforms a MFP model tailored for the particular data generating distribution. But it will be interesting if the FIM can show a comparable accuracy.

Re-iterating my first review, the study is not perfect and still exploratory. But it is interesting and certainly has a merit. I will update my score from 5 to 6.
Let's now finish with some remarks. 

- Let us note that our methodology is suitable for any choice of the coarse graining procedure. The latter is left to the practitioner, for different regions of the energy landscape can be better suited for the state definition, depending on the application at hand.
- The general applicability of our pretrained FIM was demonstrated by inferring MJPs from synthetic, simulated and experimental datasets of very different nature. However, albeit general enough to cover these cases, our **proposed training distribution should be understood as a first example**.  Our main goal was to demonstrate that designing broad synthetic training distributions can be used to train neural network models to perform (amortized) zero-shot inference of MJP in different situations. 
- Future research can involve, for example, extending our dataset to include power-law distributed rates into our datasets, and hence extend the applicability of the newly trained FIM to the case discussed above. One can also extend the pretraining dataset to include MJP of larger space state sizes, featuring connecting networks with different distributions.
As we understand it, an MJP description of glassy dynamics, as done by e.g. Bouchaud (1992), identifies the metastable states of the disordered system in question (i.e. the set of energy minima characterizing the energy landscape) with the states of the MJP. Within this description, the MJP dynamics depends on the shape of the energy landscape, inasmuch as the transition rates between states $i$ and $j$ ($f_{ij}$ in our notation) are characterized by the depth of the *trap* (that is, the height of the barrier between them).

In equations, practitioners typically write  $f_{ij} = \exp(-E_j/T)$, where $E_j$ is the trap depth and $T$ is the temperature.

What matters when applying FIM to such systems is that 
1. the number of metastable states in the system **is smaller or equal than** $C$ (the size of the largest state space in our dataset, which we set to 6 in our experiments); and that 
2. the distribution of energy barriers in the system is such that the corresponding transition rate distribution **is within our training dataset**.

Just to give an example, we can assume an exponential energy barrier distribution and find, by using the expression for $f_{ij}$ above, a rate distribution of the form

$p(f) = A  T  f^{T-1}$, for some constant $A$,

which is a power law. However, our prior rate distribution $p_{FIM}(f)$ is a (mixture of) Beta distribution(s). See e.g. line 150 or lines 473-475. Therefore, we expect FIM *trained on our current dataset*, with Beta priors, to perform poorly for glassy systems whose energy landscapes feature exponentially distributed trap depths.

We truly hope that the above answers the reviewer's question.
Thank you for the in depth responses to my questions!  I do in fact believe this extended discussion over the zero-shot terminology would go a long ways towards better communicating this concept to the reader. 

I maintain my original score.
I would like to thank the authors for their detailed explanations. My question Q3 was more about the complexity of the free energy landscapes of inference problems that can be addressed. I apologise for my abuse of language. Is the method still reliable when it comes to glassy dynamics or rugged landscapes?
*Links to references*:
- Mardt et al.(2017): https://www.nature.com/articles/s41467-017-02388-1
- Koehns et al. (2021): https://proceedings.neurips.cc/paper/2021/hash/abec16f483abb4f1810ca029aadf8446-Abstract.html
- Seifner & Sanchez (2023): https://proceedings.mlr.press/v202/seifner23a.html
- Gazzarrini et al. (2006): https://www.pnas.org/doi/full/10.1073/pnas.0600848103