category: literaturenote citekey: pessoasimulationtrainedconditionalnormalizing2025 title: "Simulation-trained conditional normalizing flows for likelihood approximation: a case study in stress regulation kinetics in yeast" authors: "Pessoa, Pedro; Martinez, Juan Andres; Vandenbroucke, Vincent; Delvigne, Frank; Pressé, Steve" year: 2025 date: 2025-00-00 2025 url: "https://arxiv.org/abs/2506.09374v1" zotero_key: Z7AJSHRG zotero_storage: NU9M4C3N collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Pessoa, Pedro"

status: converted

Simulation-trained conditional normalizing flows for likelihood approximation: a case study in stress regulation kinetics in yeast

Pedro Pessoa¹,² , Juan Andres Martinez³ , Vincent Vandenbroucke³ , Frank Delvigne³ , Steve Press´e¹,2,⁴

¹Center for Biological Physics, Arizona State University, Tempe, AZ, USA ²Department of Physics, Arizona State University, Tempe, AZ, USA ³Terra Research and Teaching Centre,

Microbial Processes and Interactions (MiPI),

Gembloux Agro-Bio Tech, University of Li`ege, Gembloux, Belgium ⁴School of Molecular Sciences, Arizona State University, Tempe, AZ, USA

Abstract

Physics-inspired inference often hinges on the ability to construct a likelihood, or the probability of observing a sequence of data given a model. These likelihoods can be directly maximized for parameter estimation, incorporated into Bayesian frameworks, or even used as loss functions in neural networks. Yet, many models, despite being conceptually simple, lack tractable likelihoods. A notable example arises in estimating protein production from snapshot measurements of actively dividing cells. Here, the challenge stems from cell divisions occurring at non-Exponentially distributed intervals with each division stochastically partitioning protein content between daughter cells, making protein counts in any given cell a function of its full division history. Such history dependence precludes a straightforward likelihood based on a (standard Markovian) master equation. Instead, we employ conditional normalizing flows (a class of neural network models designed to learn probability distributions) to approximate otherwise intractable likelihoods from simulated data. As a case study, we examine activation of the glc3 gene in yeast involved in glycogen synthesis and expressed under nutrient-limiting conditions. We monitor this activity using snapshot fluorescence measurements via flow cytometry, where GFP expression reflects glc3 promoter activity. A na¨ıve analysis of flow cytometry data ignoring cell division suggests many cells are active with low expression. However, fluorescent proteins persist and can be inherited, so cells may appear active from retaining ancestral fluorescence. Explicitly accounting for the (non-Markovian) effects of cell division reveals glc3 is mostly inactive under stress, showing that while cells occasionally activate it, expression is brief and transient.

1 Introduction

In the natural sciences, we learn to simulate dynamical models, whether deterministic or stochastic, in order to replicate experimental output and eventually predict phenomena. These models typically depend on parameters (e.g., forces or velocity fields) that can be tuned to approximate the experimental output we wish to simulate. Simulations like these have provided significant insights across scientific areas: from reactions in solution [1, 2], to gene expression [3, 4], neuronal activity [5, 6], dynamics of infection [7, 8], and even toward explaining processes involved in scientific collaboration and publishing [9, 10]. Learning parameters from data, usually termed "inverse modeling", however, demands more than the ability to simulate data from given parameters. To rigorously learn parameters, one typically formulates likelihoods — probability densities for observing data sequences given parameter sets (θ), p(data|θ). Using such likelihoods, parameters are then estimated: either by finding those parameters making the sequence of the data's likelihood largest (from which derives the notion of maximum likelihood); or by constructing posteriors over the parameters informed by the data from the likelihood [9, 11–14].

Yet, deriving likelihoods analytically is often significantly more difficult than running forward simulations. For example, while simulating chemical reactions can be straightforward via Gillespie's algorithm [1,13,15], writing down the associated likelihood involves solving the chemical master equation [4, 16, 17], requiring enumerating all possible system states feasibly visited. As the number of molecular species grows, the state space increases combinatorially, making the likelihood intractable due to large matrix exponentials [16–18]. Even this is an optimistic scenario: it applies to cases for which likelihoods can be formulated albeit at high computational cost.

Consider, for sake of example, an even more challenging case: we monitor cells at various points in time (snapshot data). From this snapshot data, our goal is to estimate the production of a particular (labeled) protein species through fluorescence. It is common to model the dynamics of this protein's number using Markovian dynamics [19–26]. For this iteration of the problem, inferring parameters is a textbook problem [13]. However, when experiments span multiple cell division cycles, the Markovian assumption fails, as many proteins have half-lives longer than a single division [27–29]. As such, reduction in protein numbers occurs, mainly, via cell division rather than degradation. Since, cell division times are not exponentially distributed [30, 31], the cell cycle cannot be modeled using Markovian dynamics. While it is possible to simulate protein levels under cell division [32–34], to our knowledge, no established inference framework allows us to write down a general likelihood for this problem.

One idea to bypass likelihood calculation is to use (forward) simulations to approximate the likelihood, as the simulated results can be understood as samples from p(data|θ), leading to what is known as simulation-based inference (SBI) [35–37]. An early implementation of SBI is approximate Bayesian computation (ABC), where the likelihood is approximated by the fraction of simulated data points falling within a threshold distance, ϵ, of observed data [38]. However, ABC presents limitations: (1) it requires new simulations for each parameter values θ, making likelihood evaluation computationally intensive; (2) it depends on arbitrary distance metrics; and (3) likelihood estimates can vary unpredictably, particularly for rare events or large datasets.

Recent advances in neural networks have inspired new approaches for SBI. Some of these consist of finding a neural network approximation of the Bayesian posterior directly, bypassing the likelihood altogether [39–42]. However, without an explicit likelihood function, we lose a key tool in validating the inferred posterior. When a likelihood is available, we can compare the approximated likelihood with histograms obtained from simulations in regions of high posterior. In some particular benchmarks [43, 44], it has been argued that posterior approximations can produce overconfident posteriors underestimating uncertainty. Without an estimated likelihood to compare against, problems like these can go unnoticed. As a result, purely posterior-focused approaches have faced criticism for limited generalizability and the lack of reliable validation strategies [43–45].

In contrast, a likelihood-focused approach provides both robustness and interpretability. For each model studied, we iteratively sample parameter values, simulate corresponding datasets, and train a neural network to map parameters to the likelihood of the data. This cycle of sampling, simulation, and refinement continues until the network accurately captures the likelihood function across the parameter space (see Fig. 1a).

Multiple neural network–based approaches have been proposed that output probability distributions and can be used to approximate likelihoods. For example, mixture density networks output the parameters of a predefined mixture distribution conditioned on model parameters [35, 36]. Energy-based models define an unnormalized probability landscape via a learned energy function and typically rely on sampling techniques to estimate likelihoods [46,47]. More recently, score-based diffusion models have been introduced, which approximate distributions by progressively corrupting data with noise and learning to reverse this process through denoising dynamics [48,49]. Other approaches shift focus from the distribution itself to learning the underlying time evolution equations, such as the chemical master equation, using neural networks thereby replacing traditional solvers during inference [50–52].

Here, we focus on normalizing flows due to their natural probabilistic interpretation, expressiveness, and ease of use for lower-dimensional applications considered herein. Normalizing flows operate by transforming a simple base distribution, typically Gaussian, through a sequence of invertible and differentiable mappings, producing a flexible and expressive target distribution. These transformations are parameterized by neural networks and trained directly on data samples to approximate unknown densities. Conditional normalizing flows [42, 53–56] extend this framework to learn probability distributions conditioned on external variables (in our case, the model parameters θ). For likelihood approximation, this means training the flow to approximate the conditional density p(data|θ) across the parameter space, providing an efficient and tractable surrogate for the true likelihood.

While conditional flows have been proposed for general density estimation tasks [42, 54], their application as intractable likelihood approximators remains relatively underexplored in scientific modeling with recent work concentrated on classification problems [55]. In this work, we broaden their scope by applying them to more general and interpretable scientific models and, in doing so, uncover the kinetics of state switching in yeast (S. cerevisiae) by accurately taking into account proteins counts across cell division.

To demonstrate the flexibility of our approach, we apply it to a hierarchy of three increasingly realistic biological models. We begin with a simple system in which a cell produces protein of interest at a constant rate and divides deterministically after fixed time intervals which we will later define as Model 1. Next, in Model 2, we extend Model 1 to include stochastic cell division times. Finally, in Model 3 we consider a realistic experimental setting where levels of activation of the glc3 gene in yeast, involved in glycogen synthesis, are monitored under nutrient-limiting conditions via fluorescent flow cytometry. By considering cell division in our inference on glc3 activity, we uncover that even under high stress conditions, glc3 remains inactive in nearly all cells with rapid switches to an upregulated state. For concreteness, a high level SBI workflow and description of our models are summarized in Fig. 1.

Figure 1: Summary of the simulation based inference with likelihood approximation for learning protein count dynamics through flow cytometry. a) Our approach involves developing a system continuously sampling parameters and simulating the dynamics, encompassing the processes underlying both the system's behavior and the measurement process, to predict the data expected from these parameters. At each new sampling step, we train a neural network to estimate the conditional probability of the observed data given the sampled parameters. b) We aim to study specific systems for application, culminating in a simulator that integrates three processes: 1) protein production, including transitions driven by phenotypical state switches [57]; 2) the splitting of proteins during cell division, a non-Markovian process that renders conventional techniques like solving the chemical master equation insufficient; and 3) the measurement process where, for this example, flow cytometry is employed to indirectly estimate the amount of a specific fluorescent protein per cell by measuring the total fluorescence intensity.

2 Methods

2.1 Simulators, latent variables and likelihoods

Simulations play a central role across scientific disciplines by enabling researchers to explore how different mechanisms influence observable outcomes. More precisely, a simulator is a computational tool that takes model parameters, denoted by θ, and uses random number generators to sample realization-specific variables, which we collectively denote by z. These variables, in turn, produce observed data points, denoted by y, that are correlated to the underlying data-generating process — such as indirect measuments of a subset of z.

A classic example of a simulator is the Gillespie algorithm [1, 13], which simulates the stochastic dynamics of chemical reaction networks. In this context, the model parameters θ represent reaction rates, while z corresponds to the full trajectory of reaction events and waiting times generated by the simulator, given these rates. In snapshot experiments, where only the state of the system at a specific time point is measured, most of the simulated trajectory z remains unobserved. We therefore refer to z collectively as latent variables. In such cases, z encodes the complete underlying stochastic history, while the observed data y captures only the state of the system at the measurement time point.

In probabilistic terms, the relationship between model parameters and observations is described by the likelihood p(y|θ). Since simulation involves latent variables, a simulator effectively provides access to samples from the joint distribution p(y, z|θ). The likelihood of the observed data can then be expressed by marginalizing over z

$$p(y|\theta) = \int dz, p(y,z|\theta) = \int dz, p(y|z,\theta), p(z|\theta). \tag{1}$$

Even when closed-form expressions for both p(y|z, θ) and p(z|θ) are available, evaluating the integral in (1) is often analytically intractable. For example, in chemical kinetics models commonly simulated with the Gillespie algorithm, the marginal likelihood p(y|θ) could, in principle, be obtained by solving the chemical master equation (CME). However, solving the CME is typically computationally intensive [16, 17, 36], and formulating the CME is only feasible for systems with Markovian dynamics. Moreover, when the latent space z is highdimensional, even numerical integration becomes prohibitively expensive.

Nevertheless, even when the marginal likelihood in (1) cannot be evaluated analytically, the simulator still allows us to generate samples from it. Provided that the simulator faithfully represents the data-generating process, drawing samples from the joint distribution p(y, z|θ) ensures that the resulting observations y follow the correct marginal distribution p(y|θ), simply by ignoring z.

This ability to sample from p(y|θ) directly enables us to approximate the likelihood without requiring explicit evaluation. In practice, we can use these samples to train a neural network to approximate the distribution p(y|θ). In the following section, we describe how this can be achieved using normalizing flows.

2.2 Normalizing flows for likelihood approximation

While simulators provide a means to generate synthetic data for known parameters, they do not, by themselves, enable inference. The key scientific challenge is to invert the simulation process: to recover the underlying parameters θ from observed data y. This requires evaluating the likelihood p(y|θ) across a range of candidate parameters. When it is infeasible to perform multiple likelihood calculations, we are inclined to use neural networks to approximate the likelihood function, treating simulations as the effective training data: by generating many (θ, y) pairs, simulators allow us to learn an approximate likelihood directly from samples in otheer words, perform SBI with neural networks [35–37].

A particularly powerful architecture for this task is the normalizing flow. Normalizing flows [47, 58–63] are a powerful tool, providing a mechanism to evaluate likelihoods by constructing a flexible family probability distributions. They work by learning differential and bijective transformations mapping a simple base distribution into the distribution we aim to approximate.

Specifically, a normalizing flow learns a transformation fϕ, parameterized by the internal parameters of the neural network, collectively denoted as ϕ, such that y = fϕ(ξ). This transformation allows us to express the probability density of y as

$$p{\rm NF}(y|\phi) = p{\xi}(\xi) \left| \det \frac{\partial f_{\phi}^{-1}}{\partial y} \right|^{-1},$$ (2)

where ξ = f −1 ϕ (y) is the inverse transformation, det ∂f ⁻¹ ϕ /∂y is the determinant of the Jacobian, and pξ(ξ) is the base distribution, usually a standard Gaussian.

For likelihood approximation, the transformation must additionally be conditioned on the model parameters θ, yielding

$$p{\rm NF}(y|\theta,\phi) = p{\xi}(\xi) \left| \det \frac{\partial f_{\phi}^{-1}(y;\theta)}{\partial y} \right|^{-1}.$$ (3)

This framework, known as conditional normalizing flows [42, 53–55], enables the model to train the neural network to learn the conditional distribution p(y|θ). In other words, to find the parameter values ϕ ∗ such that

$$p_{\rm NF}(y|\theta,\phi^*) \approx p(y|\theta).$$ (4)

To train the normalizing flow to satisfy (4), we minimize the discrepancy between the true likelihood p(y|θ) and its approximation pNF(y|θ, ϕ). A standard way to measure this discrepancy is the Kullback-Leibler (KL) divergence [64–67] defined as

$$\mathrm{KL}(\theta, \phi) = \int dy \, p(y|\theta) \log \frac{p(y|\theta)}{p_{\mathrm{NF}}(y|\theta, \phi)}. \tag{5}$$

To learn the likelihood across a range of parameter values θ, we extend the KL divergence by taking the expected value over θ from a prior distribution p(θ)

$$KL(\phi) = \int d\theta \ p(\theta)KL(\theta, \phi) = \int d\theta \ p(\theta) \int dy \ p(y|\theta) \log \frac{p(y|\theta)}{p_{NF}(y|\theta, \phi)}. \tag{6}$$

However, evaluating this expression requires explicit knowledge of the true likelihood p(y|θ), which is unavailable. Instead, we can draw a sequence of parameter samples {θs} = {θ1, θ2, . . . , θS} from the an a priori distribution, θ^s ∼ p(θ) which will later be used as the Bayesian prior, and generate corresponding synthetic data points {ys} = {y1, y2, . . . , yS} from their respective parameters, y^s ∼ p(y|θs), using the simulator, as depicted in Fig.1a, this procedure generates the set of simulated pairs (ys, θs). These samples allow one to approximate the KL divergence in (6) as a Monte Carlo integration yielding

$$KL(\theta, \phi) \approx \frac{1}{S} \sum_{s=1}^{S} \left( \log \frac{p(y_s | \thetas)}{p{NF}(y_s | \theta_s, \phi)} \right)$$

$$= C - \frac{1}{S} \sum{s=1}^{S} \log p{NF}(y_s | \theta_s, \phi).$$ (7)

where C is a constant absorbing the intractable term involving the true likelihood p(ys|θs), which does not depend on the neural network parameters ϕ and therefore does not affect optimization.

As a result, minimizing the KL divergence is equivalent to minimizing the empirical negative log-likelihood of the simulated data

$$\mathcal{L}(\phi) = -\frac{1}{S} \sum{s=1}^{S} \log p{\text{NF}}(y_s | \theta_s, \phi)$$ (8)

and obtaining

$$\phi^* = \arg\min_{\phi} \mathcal{L}(\phi). \tag{9}$$

This approach allows us to sidestep the need to compute the true KL divergence, relying instead on evaluating the conditional normalizing flow probability pNF(ys|θs) for simulated pairs. As training depends on these simulated (θs, ys) pairs, we adopt a dynamic resampling strategy during training to avoid overfitting to a finite number of simulations. As illustrated in Fig.1a, each training cycle consists of sampling a new set of parameter values {θs} from the prior and running the simulator to generate corresponding synthetic data {ys}. These newly generated pairs are then used to calculate the loss function defined in (8), treating (θs, ys) as training examples. This cycle (sampling parameters, simulating observations, and updating ϕ in order to minimize (8)) is repeated throughout the training to ensure that the neural network learns a smooth and generalizable approximation of the likelihood across the full parameter space. The specific details of the resampling frequency, dataset sizes, and replacement strategies depend on the model under study and are described in Sec.3.

2.3 Bayes' theorem and parameter learning

Having trained the normalizing flow ^pNF(y|θ, ϕ^∗ ) to approximate the likelihood p(y|θ), we now turn to the problem of inferring model parameters θ from real observed data. We denote the observed dataset as {y ⁿ} ⁼ {^y 1 , y² , . . . , y^N }, where each ^y n represents an experimental measurement. This observed data should not be confused with the simulated samples used during training, denoted as ys.

From Bayes' theorem, the distribution for θ given {y ⁿ}, the posterior, is

$$p(\theta|{y^n}) = \frac{p({y^n}|\theta)p(\theta)}{p({y^n})} \propto p({y^n}|\theta)p(\theta).$$ (10)

where p(θ) is the prior over the parameters, and p({y ⁿ}) is the evidence, independent of ^θ, thus allowing us to write the posterior up to proportionality as above.

Assuming that the data points are conditionally independent given θ, the total likelihood factorizes as p({y ⁿ}|θ) = ^Q^N ⁿ=1 p(y n |θ). We can now replace each likelihood factor with the output of the trained normalizing flow, setting p(y n |θ) ≈ pNF(y n ^|θ, ϕ^∗ ). Substituting this into the factorized likelihood, the posterior becomes

$$p(\theta|{y^n}) \propto \left(\prod{n=1}^N p{\rm NF}(y^n|\theta,\phi^*)\right) p(\theta),$$ (11)

In the following sections, we apply this framework to specific models of interest. For each case, we define the forward process that generates synthetic data, train a conditional normalizing flow to approximate the likelihood, and perform Bayesian inference based on the observed experimental data.

3 Results

Here, we evaluate the performance of our likelihood approximation framework, in a sequence of models of increasing biological complexity:

Model 1: Regular cell division (Sec. 3.1) This model admits an exact analytical expression for the steady-state likelihood serving as a reference case, and allowing us to quantitatively compare our neural network-based approximation to the exact solution.
Model 2: Stochastic cell division (Sec. 3.2) We introduce stochasticity in division times rendering the likelihood analytically intractable. We validate the neural network approximation by comparing the predicted likelihood at given parameters to simulation results.
Model 3: Stochastic cell division with two-state regulation and uneven partitioning (Sec. 3.3) This model introduces additional biological complexity, including state switching regulation, asymmetric division, and indirect fluorescence measurements. It is fully simulation-based and designed to reflect realistic S. cerevisiae protein dynamics, ultimately capable of replicating the experimental conditions recovered in [57] for which we present the inference results in Sec. 3.4.

3.1 Model 1 — Regular cell division model

The first model describes a stochastic system in which protein production occurs at a constant rate β, and cells divide at regular intervals T. Importantly, this model admits an exact analytical expression for the steady-state likelihood of protein counts, allowing us to directly compare the neural network-based likelihood approximation to the exact result.

The system includes the following sources of stochasticity: (1) Protein synthesis is modeled as a birth process, meaning the number of proteins produced within a time interval ∆t follows a Poisson distribution with mean β∆t; (2) Cell division events are asynchronous, meaning the cell can be found in any stage of the cell cycle, independent of other cells in the environment. We denote by τ ∈ (0, T) as the time since the last division event. Each cell has a τ uniformly distributed at any point in the cell cycle, p(τ ) = 1/T for τ ∈ (0, T); (3) The partitioning of proteins between the two daughter cells is idealized as Binomial for this simple model. Put differently, each protein is transferred to one of the daughter cells with probability 1/2. Thus, if the parent cell contains G^T proteins preceding division, the number of proteins allocated to one daughter cell follows a Binomial distribution with G^T trials and a success probability of 1/2.

3.1.1 Steady-state distribution

Here, we show that under the sources of stochasticity enumerated above (points 1–3), it is possible to derive a steady-state distribution for the number of proteins, denoted by G. The stochastic dynamics corresponding to conditions 1 and 3 can be formalized mathematically as

$$G - G_0|\beta, \tau \sim \text{Poisson}(\beta\tau),$$ (12a)

$$G_0|\beta \sim \text{Binomial}(G_T, 1/2),$$ (12b)

where ∼ denotes "sampled from" or "distributed as". Above, G⁰ and G^T refer to the protein counts immediately after and immediately before a division event, respectively. Since we later marginalize over τ , we simplify the notation by treating G as the protein count at an arbitrary stage in the cell cycle, with τ being a latent variable.

A distribution consistent with the conditions in (12) is

$$G|\beta, \tau \sim \text{Poisson}(\beta(1+\tau)).$$ (13)

From the expression above, we can recover the distributions of G⁰ and G^T by evaluating G^τ for any value of τ . In particular, setting τ = 0 and τ = 1 yields:

$$G_0|\beta \sim \text{Poisson}(\beta),$$ (14a)

$$G_T|\beta \sim \text{Poisson}(2\beta).$$ (14b)

As a sanity check, these results are consistent with the constraints in (12). For instance, as the sum of two independent Poisson random variables is itself Poisson distributed, with the mean equal to the sum of the individual means, the combination of (13) and (14a) satisfies (12a). Second, a well-known property of Poisson distributions is that partitioning a Poissondistributed variable via a Binomial process results in another Poisson-distributed variable, where the mean is scaled by the Binomial success probability. Thus, (14a) and (14b) are consistent with (12b).

Now, in order to have a closed-form expression for p(G|β), we observe that the partitioning of proteins follows a Binomial between daughter cells), we observe cells at uniformly distributed times within the cell cycle, allowing us to marginalize over τ in (13) as

$$p(G|\beta) = \int d\tau \, p(\tau) p(G|\beta, \tau)$$

$$= \int_0^T d\tau \, \frac{1}{T} \frac{e^{-\beta(1+\tau)} (\beta(1+\tau))^G}{G!}$$

$$= \frac{\gamma(G+1, 2\beta) - \gamma(G+1, \beta)}{\beta T}.$$ (15)

where γ(a, x) represents the incomplete gamma function

$$\gamma(a,x) = \frac{1}{\Gamma(a)} \int_0^x dt \ t^{a-1} e^{-t}.$$ (16)

This analytical steady-state distribution was obtained from a model where proteins disappear through partitioning at cell division rather than through degradation. This contrasts with the classical birth–death framework commonly used in the literature [19–26], where individual proteins degrade at a constant rate γ. Supposing a constant protein production rate β, the steady-state protein distribution in the birth–death model is Poisson with mean β/γ. To highlight the dramatic qualitative differences between birth-death process versus our cell-division process, Fig. 2 shows a comparison between simulations of a birth–death process with parameters β = 150γ, and the cell-division process described here with β = 100/T, adjusted to achieve the same expected steady-state mean according to (15).

Having established the steady-state distribution for protein counts in (15), along with a qualitative comparison to the classical birth-and-death framework illustrated in Fig. 2, we have a basis to which to compare the results of our likelihood approximation scheme for this model.

3.1.2 Training the neural network

Our goal is to train a neural network to approximate the likelihood, i.e., the steady-state distribution of protein counts. Having an exact analytical solution for the steady-state distribution in (15), we use the regular cell division process model as an initial test case to train

Figure 2: Why the Markovian approach is insufficient to learn protein production counts across cell division. RNA production is often modeled as a birth-death process 4, 17, where constant production leads to a Poisson steady-state distribution (top right). While previous literature extend these Markovian models to account for protein production [19–23], marker proteins have lifetimes on the order of hundreds of hours, thus making normal elimination by mean of degradation negligible as compared to the reduction caused by cell division. In the case where protein production occurs at a constant rate and cell divisions happen deterministically at equal intervals (described in Sec. 3.1 and simulated in the bottom left), the steady-state distribution becomes an incomplete Gamma distribution (15) (bottom right), which is qualitatively different from the Poisson distribution recovered for birthdeath process. However, more complex scenarios, such as variability in division times (introduced in Sec. 3.2) or non-constant production rates (discussed in Sec. 3.3) complicate the model even further, necessitating a simulation-based approach to accurately model protein counts.

and evaluate the network, before moving on to models where the exact distribution is not available.

When training the network, we work in log space for both the parameters and the protein counts G. Specifically, we sample values of ln β, transform back to β, run the simulator with β to obtain protein counts G, and then take the logarithm again to obtain ln G. Thus, ln G takes the role of the observations (y in Sec. 2), while ln β takes the role of the parameters (θ in Sec. 2). Training the network consists of finding ϕ ∗ such that

$$p(\ln G|\ln \beta) \approx p_{\rm NF}(\ln G|\phi^*, \ln \beta).$$ (17)

If the original distribution over G is needed, it can be recovered as

$$p(G|\beta) = p(\ln G|\ln \beta) \frac{d \ln G}{dG}$$

$$\approx \frac{1}{G} p_{NF}(\ln G|\phi^*, \ln \beta).$$ (18)

Having outlined the mapping between simulation outputs and network inputs in log space, we now describe how the training data is generated. The training process begins by sampling ln β, with β expressed in units of 1/T, from a Gaussian distribution with mean 7 and variance 1. This prior choice is only used for this example. When training a network to deal with experimental data, we train a new network using data informed parameter estimates. For each sampled ln β, we compute the corresponding β and simulate the dynamics of protein production and division. This consists of sampling an initial latent variable, assuming that the cell's last division occurred at a random time drawn uniformly between 0 and T. Each simulated cell is then evolved for a duration equivalent to 100 cell division intervals, ensuring convergence to the steady-state protein distribution. Between divisions, the number of proteins generated is sampled from a Poisson distribution according to (12a) and at the end they are divided as a Binomial sample as in (12b). Further details on how we efficiently parallelize multiple simulations on GPUs are available in our GitHub repository [68].

The training dataset consists of about a million, 220, pairs of (ln β, ln G) values, divided into 16 batches. For each batch, we compute the loss function as defined in (8) and perform one optimizer step to update the network parameters. A full pass over all batches, going once through the entire dataset, defines one epoch. The network is trained over many epochs, progressively minimizing the loss and resulting in optimized parameters ϕ ∗ , as described in Sec. 2.2.

To prevent overfitting, at each epoch we regenerate 2¹⁴ new pairs to replace part of the dataset. This practice allows the network to be exposed to different samples in each training iteration, thereby promoting better generalization and reducing the risk of memorizing specific data points. Additionally, we shuffle the dataset indices after every epoch to create varied batch compositions and to further expose the network to a broad range of parameter values during training.

We made the code used to generate the dataset, train the network, and perform all subsequent analyses publicly available on our GitHub repository [68]. Our implementation builds on the normflows package [62] for normalizing flows. Training was performed on a machine equipped with an 8-core Intel Core i7-7700K CPU, 64 GB of system RAM, and an NVIDIA GeForce RTX 3070 GPU with 8 GB of dedicated memory. A typical training run, consisting of 1000 epochs, takes approximately 1.5 hours. These hardware specifications apply to all other benchmarks presented in this work.

3.1.3 Inference

We present the results obtained using a normalizing flow trained over 2000 epochs in Fig. 3. Fig. 3a summarizes the model. In Fig. 3b, we compare the neural network's approximation of the likelihood function for selected β values to both the exact solution and data simulated using the corresponding β values (orange histogram). Notably, the neural network was trained on a diverse range of β values sampled across the parameter space, rather than specifically on the selected values shown. This comparison demonstrates the network's ability to accurately approximate the likelihood at a given point in the parameter space, even without extensive training data at that specific β value.

Finally, in Fig. 3c, we illustrate how the neural network–based likelihood approximation can be used for Bayesian inference to estimate β from synthetic data. We consider sets of N independent snapshots observations Gⁿ collectively denoted simply as{Gn} ⁼ {G¹ , G² , . . . , G^N }, and compute the posterior distribution over a grid of ^β values, with the total likelihood approximated as

$$p({G^n}|\beta) \approx \prod{n=1}^{N} \frac{1}{G^n} p{NF}(\ln G^n | \phi^*, \ln \beta),$$ (19)

where pNF is the output of the trained normalizing flow model with parameters ϕ ∗ , obtained from the training scheme described in Sec. 3.1.2, and capture the conditional distribution of ln G given ln β. The posterior is then given by

$$p(\beta|{G^n}) \propto p({G^n}|\beta) p(\beta), \tag{20}$$

Figure 3: Results for Model 1. a) Schematic of the model described in Sec. 3.1. The model involves a single parameter, β, representing the protein production rate. Proteins are diluted during cell division, which occurs at fixed intervals of duration T. b) Neural network approximation of the likelihood (blue) for selected β values, compared to data simulated using these values as ground truth (orange histogram) and the exact likelihood (15) (black). The neural network accurately captures the likelihood across a broad range of β values. c) Application of Bayesian inference to estimate β from synthetic data. The posterior distribution tightens as the dataset size increases, reflecting improved parameter estimation with more observations. Even with a modest number of observed cells, the inference achieves reasonable accuracy in this example.

and the proportionality constant is obtained via numerical integration over the grid of β values. This posterior reflects how uncertainty about β narrows as more data is incorporated. In Fig. 3c, we show posterior distributions for increasing dataset sizes N, highlighting the gain in precision with larger sample sizes.

3.2 Model 2 — Stochastic cell division

In the previous section, we considered a system in which protein production occurs at a constant rate β, and cell division takes place at regular time intervals T. We now extend this model to account for the fact that division timing is inherently stochastic, due to both environmental and intracellular fluctuations. In many practical scenario, such as chemostat experiments, the average division time, T, can be independently estimated or controlled, as further discussed in Supplemental Information S1. We therefore treat T as known and serving as a natural unit of time.

While a comprehensive treatment of mechanisms causing variability in cell division timing is beyond the scope of the present work [30, 31, 69], we adopt a phenomenological approach in which division times are modeled as the sum of a large number of sequential, memoryless steps [70–72]. This framework leads naturally to the use of Gamma distributions.

Thus, we model division times as Gamma distributed in units of the average division time, T, that is

$$\frac{T}{\overline{T}} \sim \text{Gamma}\left(\frac{1}{\sigma^2}, \sigma^2\right),$$ (21)

where we introduced a new parameter, σ, equivalent to the standard deviation of T /T. This formulation allows us to flexibly capture biologically realistic fluctuations in division timing while maintaining analytical control over the statistical properties of the system.

In this framework, the parameters to be inferred include the protein production rate, β, and the standard deviation of the cell division time, σ, which captures the stochasticity of the cell cycle as defined in (21). As in Model 1, the number of proteins produced during a time interval ∆t between cell divisions is Poisson distributed with mean β∆t and upon each division event the number of proteins in the cell is reduced according to a Binomial distribution with success probability 1/2. A summary of the model is presented in Fig. 4a.

For this model, deriving an exact analytical expression for the steady-state protein distribution, analogous to (15), becomes infeasible. Consequently, we rely on simulation-based methods combined with normalizing flow-based likelihood approximation to perform parameter inference, following the framework introduced in the previous section.

3.2.1 Training the neural network

The training process follows a structure similar to that used in the deterministic case of Model 1 (Sec. 3.1.2), but now incorporates stochastic division times. Specifically, the neural network is trained to approximate the probability distribution p(G|β, σ). During training, each parameter set, θ = (ln β, ln σ), is sampled, where ln σ is drawn from a Gaussian distribution with mean −2.3 and standard deviation 0.5, and ln β is drawn from a Gaussian with mean 7 and standard deviation 1 (these priors are used only for training on synthetic data). As before, the inputs to the network are ln G and θ, and the output is the learned likelihood function pNF(ln G|ϕ ∗ , θ).

As in the deterministic case in Model 1, we again generate a dataset of 2²⁰ pairs of (θ, ln G), split into 16 batches. The network is trained iteratively using the negative log-likelihood loss function, as defined in (8). At each epoch, 2¹⁴ new pairs are sampled to replace a fraction of the dataset, and the data indices are shuffled to ensure diverse batch compositions. In the

Figure 4: Results for Model 2. a) Schematic of the model described in Sec. 3.2. We learn the pair of parameters θ = {β, σ} representing the protein production rate and cell division time standard deviation. b) Neural network approximation of the likelihood (blue) for selected θ values, compared to data simulated using these values as ground truth (orange histogram). The neural network accurately captures the likelihood across a broad range of β values. c) Application of Bayesian inference to estimate both parameters from synthetic data. The posterior distribution tightens as the dataset size increases, from 512 to 2048 cells, reflecting improved parameter estimation with more observations.

same machine specifications as described in Sec. 3.1.3, the training consisting of 1000 epochs using the code available on GitHub took approximately 2 hours.

3.2.2 Inference

We now evaluate the performance of the trained normalizing flow model in approximating the likelihood function in Model 2, modeling stochastic division dynamics. The results are summarized in Fig. 4.

Fig. 4a provides an overview of the model described in Sec. 3.2, highlighting the two parameters to be learned (β and σ). In Fig. 4b, we compare the neural network's likelihood approximation and simulations at these parameter values. Finally, in Fig. 4c, we apply Bayesian inference to estimate both β and σ from synthetic data. The posterior distributions illustrate how the inferred values become more constrained as more observations are incorporated. This result highlights the effectiveness of the neural network-based likelihood approximation in supporting inference for stochastic division models.

3.3 Model 3 — Stochastic cell division with two-state regulation and uneven partitioning

We now extend our model by introducing additional layers of biological complexity. With the intent of applying this method to flow cytometry data of S. cerevisiae, we incorporate a twostate regulatory mechanism governing protein production, uneven cell division (as observed in yeast), and the indirect measurement of protein content via flow cytometry. Due to the difficulties introduced by these additional biological features, Model 3 does not admit an analytical solution for p(G|θ).

In this extended framework, cells transition between two internal states: an active gene state, where protein production occurs at rate β, and an inactive state, where the protein of interest, whose counts are observed indirectly through fluorescence, is not produced. Simulating this model requires keeping track of both the number of proteins and the latent internal gene state over time. The time until the next gene state switch is exponentially distributed, with transition rates given by the activation rate λact (inactive to active) and the inactivation rate λina (active to inactive). The parameters to be inferred now include (β, σ, λact, λina).

As with the previous model, cell division times T are sampled from (21), while the times to state switching events are sampled from Exponential distributions determined by the corresponding transition rates. We assume that the production of the protein of interest occurs only during the active gene state. Thus, within any time interval ∆t, defined as the time between cell division or state switching events, the number of proteins of interest produced is modeled as a Poisson random variable with mean β∆t.

Unlike in Model 2, we do not assume that cells divide symmetrically. Rather, we model the volume fraction ρ of the daughter cell relative to the total volume as a Beta-distributed random variable. Drawing on previous chemostat fermentation data at various dilution rates reported by some of us [57], we have ρ ∼ Beta(6, 14), corresponding to a mean of 0.3 and a standard deviation of 0.1. During each division event in the simulation, the cell is randomly assigned to be either a mother or a daughter and the number of proteins inherited is determined by a Binomial distribution with success probability ρ (if the cell is assigned as the daughter) or 1 − ρ (if assigned as the mother). In other words, proteins are partitioned independently, with inheritance probabilities proportional to the corresponding cell's volume fraction.

Although previous models assumed direct count measurements for the protein of interest, here we model protein abundance as being correlated with fluorescence intensity. In this framework, we must account for autofluorescence: a background signal emitted by cells at the same frequency as the protein marker.

The autofluorescence intensity per cell, denoted by Iaf, is approximated using a fixed normalizing flow trained on single-cell fluorescence measurements collected under minimal stress conditions (see Supplemental Information S2). This allows us to model the observed fluorescence I as the sum of the protein fluorescence signal I^g and the autofluorescence contribution:

$$I = I_{\rm af} + I_g. \tag{22}$$

To incorporate fluorescence measurements into our model, we first sample Iaf from the pre-trained normalizing flow. Next, the fluorescence signal I^g of the protein of interest is generated according to:

$$I_g \sim \text{Gaussian}(\xi G, \xi G),$$ (23)

where ξ is a calibration factor reflecting the conversion of fluorescent particles (G) into total fluorescence intensity (Ig). ξ was determined from bead experiments (see Supplemental Information S3).

Finally, the observed fluorescence I is obtained as (22), and the neural network takes as input the log-transformed intensity ln I along with the logarithm of the parameters ψθ. The output is the learned likelihood function pNF(ln I | ϕ ∗ , θ).

3.3.1 Training the neural network

The training process follows a structure similar to the previous examples in Sec. 3.1 and 3.2. Each parameter set is sampled as θ = (ln β, ln σ, ln λact, ln λina), where ln σ is drawn from a Gaussian distribution with mean −2.3 and standard deviation 0.5, ln β is sampled from a Gaussian distribution with mean 7 and standard deviation 1, ln λact with mean 1 and standard deviation 1.5, and ln λina with mean −1 and standard deviation 1.5. These priors apply only to synthetic data; later when applying to real data, we retrain the network using calibrated estimates as detailed in Sec. 3.4.

To ensure robust training, we again generate a dataset of 2²⁰ pairs of (ln I, θ) split into 16 batches. The network is trained iteratively using the negative log-likelihood loss function defined in (8). At each epoch, 2¹⁴ new pairs are sampled to replace a fraction of the dataset. Training the model for 1000 epochs using the code available on GitHub took approximately 2.5 hours.

3.3.2 Inference

We now evaluate the performance of the normalizing flow model in approximating likelihoods and enabling Bayesian inference in this more complex setting, which includes a two-state gene activation regulatory mechanism, uneven protein partitioning, and fluorescence-based measurements.

Fig. 5a provides an overview of this extended model. Unlike previous sections, where cell division was modeled as a symmetric split (bacterial division) this model accounts for the asymmetric division of yeast, which occurs through budding. In this process, daughter cells inherit different amounts of protein depending on their size at division.

The increased number of parameters (β, σ, λact, λina) makes direct posterior evaluation across a range of values in the region of interest computationally impractical. While this approach was feasible in the example of Sec. 3.1, which involved a single parameter β, and in Sec. 3.2, which considered two parameters β and σ, it becomes intractable in higher dimensions. We therefore employ a posterior sampling though an Markov chain Monte Carlo MCMC method with an adaptive proposal strategy inspired by [73], which dynamically refines

Figure 5: Results for Model 3. a) Schematic of the model described in Sec. 3.3. Including the budding-type division that is common in S. cerevisiae. We learn the set of parameters θ = {β, σ, λact, λina} representing the protein production rate, cell division time standard deviation, activation rate, and deactivation rate. b) Neural network approximation of the likelihood function for the total flow cytometry intensity I, in arbitrary florescence units, shown in blue for selected parameter values θ. The orange histograms represent samples drawn from simulations of the model using these parameter values. c) Bayesian inference results for parameter estimation. Each column corresponds to a different parameter, with inferred posterior distributions shown as blue histograms obtained via MCMC. The red dashed lines indicate the ground truth values used to generate the synthetic data. These posteriors are obtained with 2048 datapoints. Each row in c) corresponds to the same ground truth values used in the respective panel in b).

the sampling distribution for efficiency. At each MCMC step, we evaluate the posterior by substituting the true likelihood with its approximation learned via a normalizing flow, as described in Sec. 3.1.2. Further details on this implementation are provided in Supplemental Information S4.

In Fig. 5b, the learned likelihood functions exhibit qualitatively distinct patterns depending on the parameter values. Despite this variability, a single neural network successfully approximates the likelihood across a broad range of conditions, demonstrating its flexibility and accuracy.

Furthermore, Fig. 5c shows that the inferred posterior distributions are well-separated for different ground truth values, confirming the neural network's ability to resolve distinct kinetic regimes. This separation highlights the effectiveness of the method in distinguishing parameter variations, even within a complex biological system.

With this trained network incorporating realistic biological features relevant to S. cerevisiae, we now apply it to experimental flow cytometry data to infer in vivo protein kinetics.

3.4 Learning kinetics from flow cytometry of S. cerevisiae

As a real-world application of the model developed in the previous sections, we analyze flow cytometry data from a study previously conducted by some of us [57], which examined the general stress response of S. cerevisiae under continuous cultivation. This yeast species is known to employ bet-hedging strategies to balance growth and stress resistance [74]. As a result, even under optimal growth conditions, genetically identical cells may exhibit divergent phenotypes with some continuing to grow while others invest in stress resistance, thus leading to significant heterogeneity in gene expression.

To monitor this heterogeneity, a transcriptional reporter was developed using GFP under the control of the promoter of the glc3 gene, which is involved in glycogen synthesis and is known to be modulated by nutrient limitation [57, 75]. This reporter was previously shown to be a marker of bet-hedging in yeast [74].

To systematically vary nutrient availability, chemostat experiments were conducted at three dilution rates: 0.12, 0.23, and 0.33h −1 . At the lowest rate, nutrient limitation is high and the glc3 reporter is strongly expressed. At the highest rate, nutrients are abundant, suppressing reporter expression. The intermediate condition is expected to yield a heterogeneous mix of stressed and growing cells. During these cultivations, flow cytometry was used to quantify GFP expression in individual cells [57, 76].

Here flow cytometry analyzes individual events, where each event corresponds to a single particle, usually a single cell, passing through a laser beam. For this analysis, we focus on forward scatter (correlating with cell size [77]), side scatter (providing orthogonal size information), and FL1 fluorescence intensity quantifying the amount of GFP, a proxy for the expression of the glc3.

Flow cytometry data, however, can be impacted by artifacts. Debris and non-cellular particles may be misclassified as valid events (cells). To ensure data quality, we apply a preprocessing pipeline removing those artifacts. Following recomendation from Flow cytometer manuals [78], we first impose a lower cutoff on forward scatter is used to exclude non-cellular particles. Then, a support vector machine classifier [79] is employed to detect and remove outliers, as validated in other studies of both flow and mass cytometry [80, 81]. This procedure discards only a small fraction of events but significantly improves the reliability of downstream inference. Full details are provided in Supplemental Information S5.

The resulting fluorescence intensity distributions are shown in Fig. 5a. These distributions support the assumption, further discussed in Supplemental Information S6 that only one gene state is primarily responsible for protein production.

Figure 6: Inference of protein kinetics from real flow cytometry data. (a) Fluorescence intensity, in arbitrary florescence units, distributions for yeast populations at different dilution rates. Since lower dilution rates correspond to higher nutrient limitation and thus higher stress, the high stress condition (0.12 h −1 ) shows significantly elevated fluorescence compared to the low stress condition (0.23 h −1 ), while the control (0.33 h −1 ), used to calibrate autofluorescence, exhibits only background signal. (b) Posterior distributions of the inferred kinetic parameters obtained through MCMC. The model assumes a constant protein production rate (β), while gene activation (λact), inactivation (λina), and relative division-time variability (σ) are different at each stress condition. (c) Estimated fraction of time spent in the active state, computed via 24. Cells in high stress are active ≈ 5% of the time, versus ≈ 2% in low stress. Although activation is brief, proteins persist across several cell cycles, explaining the high observed fluorescence.

In the chemostat, the dilution rate determines the rate of nutrient resupply and thus serves as a proxy for cellular stress. Lower dilution rates increase nutrient limitation, activating stress responses. Among the tested conditions, 0.33 h −1 is the highest dilution rate at which a stable population could be maintained. As nutrient levels are high and stress is minimal, this condition is used to measure autofluorescence, as detailed in Supplemental Information S2. The dilution rate of 0.23h −1 represents low stress, while 0.12h −1 corresponds to high stress and prolonged nutrient limitation.

Since the cell division time in a chemostat is the inverse of the dilution rate, see Supplemental Information S1, we analyze data in rescaled time units such that the average cell cycle division time is the unit of time. Inference is performed in these units and later converted back to hours, yielding average division times of approximately 4.35h (low stress) and 8.33h (high stress).

As shown in Fig. 6a, fluorescence intensity is orders of magnitude higher in the high stress condition compared to low stress and the autofluorescence control. At first glance, this may suggest that cells under high stress remain constantly in the active state, producing fluorescent proteins continuously. However, a deeper analysis of the inferred kinetic parameters reveals a more nuanced picture.

In the network training performed on synthetic data, we sampled model parameters θ from a predefined distribution. However, without prior knowledge of which parameter values are

actually compatible with the real data, this choice may lead to suboptimal training. Ideally, we want to train the neural network around parameter regions where the likelihood is expected to be high. This approach ensures that the network sees more examples of datasets in regions similar to the real data, thereby leading to better approximations of the likelihood function. At the same time, it is important to maintain sufficient flexibility, meaning including a wide enough range of parameter values to capture possible biological variability and account for uncertainties. To achieve this balance in way compatible with the data in Fig. 6a, we train the neural network on data generated from parameters calibrated to the real experiment using a more naive form of SBI, as detailed in Supplemental Information S7.

We then perform MCMC inference, similarly to how it was done in Sec. 3.3.2 (with further details in Supplemental Information S4), and the posterior obtained is presented in Fig. 6b. The estimated protein production rate, β, is assumed unknown and to be learned but is assumed to be the same across both stress conditions. This is because β is set by the genetic machinery (e.g., mainly related to the ribosome binding) and is not directly tied to environmental factors, while the other parameters: cell division variability (σ), gene activation rate (λact), and inactivation rate (λina) are condition dependent and differ across high and low stress conditions.

From these inferred parameters, we compute the fraction of time spent in the active state, given by

$$f{\rm active} = \frac{\lambda{\rm act}}{\lambda{\rm act} + \lambda{\rm ina}}.$$ (24)

The resulting estimates, shown in Fig. 6c, reveal that cells under high stress actually spend only about 5% of the time in the active state, compared to 2% in the low stress condition. This finding suggests that under high stress, cells are not persistently active but instead switch into the active state only briefly. However, due to the slow degradation of the produced proteins, fluorescence levels remain elevated across multiple cell cycles, leading to the observed high fluorescence intensity in Fig. 6a. This underscores the importance of modeling cell division explicitly and the necessity of the current framework: a na¨ıve interpretation of raw fluorescence would incorrectly suggest that a large fraction of cells were active. In contrast, capturing cell division as part of the inference reveals that gene activation is rare and shortlived. Our framework thereby helps us disentangle inherited fluorescence from true expression dynamics.

4 Conclusion

Inference is a critical step in going from data to a parametrized model. Yet physics-inspired inference is limited by our ability to mathematically formulate likelihoods, the probability of observing the data given the model. This often reduces to solving the coinciding master equation [3, 4, 16, 17], Fokker-Planck [11, 82, 83] and other often equations which often cannot be solved or easily formulated though the process can be simulated.

In this work, we demonstrated how neural networks can be leveraged as approximators to likelihoods using data from forward simulations in inference problems relevant not only to simple to simulate yet pathologically non-Markovian dynamics but also to inference problems where the likelihood is computationally expensive. In such settings, the method presented here provides a practical surrogate preserving the statistical framework of likelihood-based inference while substantially reducing computational costs.

While we have opted for normalizing flows in this work due to their straightforward probabilistic formulation, though other neural network architectures whose outputs are tractable probability distributions (such as score-based diffusion models [48]) may also have been used. For example, diffusion models have been applied to infer Boltzmann weights for proteins [84], offering a promising alternative generative modeling framework.

As an illustrative example, we applied this approach to a non-Markovian biological system involving protein production kinetics, where protein abundance is reduced through cell division, thereby making the instantaneous protein count within a cell dependent on its complex cell-division history.

This feature essentially undermines simplifying assumptions used thus far to make this inference problem tractable such as treating cell division on protein abundance is frequently modeled as an additional chemical reaction [19–26] which assumes that reduction in protein numbers arises solely from individual degradation rather than partitioning during inheritance, in a non-Markovian cell cycle. As shown in Fig. 2, such approximations can lead to qualitatively incorrect likelihood distributions. In contrast, our method models cell division times explicitly and performs likelihood-based inference directly from forward simulations, thereby avoiding these oversimplifications.

We applied this likelihood approximation technique to both deterministic (Model 1) and stochastic (Model 2) models of cell division time, as shown in Figs. 3 and 4, respectively. In the deterministic case, the method accurately captured the steady-state protein distribution described in (15), demonstrating its ability to learn complex and otherwise intractable likelihood functions. When extended to stochastic division times, the method also proved effective in multi-parameter inference, successfully estimating both the protein production rate, β, and the division time variability, σ.

We further incorporated additional biological complexity, including two-state gene regulation and asymmetric cell division, generalizing inference to more realistic scenarios such as those found in flow cytometry datasets of S. cerevisiae. The extended model allowed us to infer phenotypic gene activation and inactivation rates, λact and λina, from synthetic data (Fig. 5).

We then applied this framework to analyze real flow cytometry data from S. cerevisiae cultured in chemostats at different dilution and thus nutrient concentration rates (Fig. 6). The results show that contrary to expectation — and although fluorescence intensity is significantly higher under high-stress conditions — cells rarely find themselves in a glc3 upregulated state [57,76] instead leveraging gene product inherited through parent cells. Interpreting these results without accounting for inheritance across the cell cycle would lead to fundamentally different conclusions, suggesting that glc3 is frequently upregulated ignoring altogether gene products inherited through division.

Acknowledgments

All authors thank Roberto Covino for inspiration in simulation-based inference and encouraging us to think deeply about this new approach. SP acknowledges support from the NIH (R35GM148237), ARO (W911NF-23-1-0304), and NSF (Grant No. 2310610).

References

[1] D. T. Gillespie, Exact stochastic simulation of coupled chemical reactions, The Journal of Physical Chemistry 81 (1977) 2340–2361.
[2] D. Schnoerr, G. Sanguinetti, R. Grima, Approximation and inference methods for stochastic biochemical kinetics—a tutorial review, Journal of Physics A: Mathematical and Theoretical 50 (2017) 093001. doi:10.1088/1751-8121/aa54d9.
[3] J. Szavits-Nossan, R. Grima, Solving stochastic gene-expression models using queueing theory: A tutorial review, Biophysical Journal 123 (2024) 1034–1057. doi:10.1016/j.bpj.2024.04.004.
[4] Z. Kilic, M. Schweiger, C. Moyer, D. Shepherd, S. Press´e, Gene expression model inference from snapshot RNA data using Bayesian non-parametrics, Nature Computational Science 3 (2023) 174–183. doi:10.1038/s43588-022-00392-0.
[5] N. Kriegeskorte, P. K. Douglas, Cognitive computational neuroscience, Nature Neuroscience 21 (2018) 1148–1160. doi:10.1038/s41593-018-0210-5.
[6] I. Blundell, R. Brette, T. A. Cleland, T. G. Close, D. Coca, A. P. Davison, S. Diaz-Pier, C. Fernandez Musoles, P. Gleeson, D. F. M. Goodman, M. Hines, M. W. Hopkins, P. Kumbhar, D. R. Lester, B. Marin, A. Morrison, E. M¨uller, T. Nowotny, A. Peyser, D. Plotnikov, P. Richmond, A. Rowley, B. Rumpe, M. Stimberg, A. B. Stokes, A. Tomkins, G. Trensch, M. Woodman, J. M. Eppler, Code generation in computational neuroscience: A review of tools and techniques, Frontiers in Neuroinformatics 12 (2018) 68. doi:10.3389/fninf.2018.00068.
[7] M. S. Jalali, C. DiGennaro, A. Guitar, K. Lew, H. Rahmandad, Evolution and reproducibility of simulation modeling in epidemiology and health policy over half a century, Epidemiologic Reviews 43 (2021) 166–175. doi:10.1093/epirev/mxab006.
[8] W. Waites, M. Cavaliere, D. Manheim, J. Panovska-Griffiths, V. Danos, Rulebased epidemic models, Journal of Theoretical Biology 530 (2021) 110851. doi:10.1016/j.jtbi.2021.110851.
[9] P. Pessoa, S. Press´e, How many submissions are needed to discover friendly suggested reviewers?, PLOS ONE 18 (2023) 1–17. doi:10.1371/journal.pone.0284212.
[10] D. Seˇselja, Agent-based models of scientific interaction, Philosophy Compass 17 (2022) ˇ e12855. doi:10.1111/phc3.12855.
[11] J. S. Bryan, I. Sgouralis, S. Press´e, Inferring effective forces for Langevin dynamics using Gaussian processes, The Journal of Chemical Physics 152 (2020) 124106. doi:10.1063/1.5144523.
[12] J. S. Bryan, P. Basak, J. Bechhoefer, S. Press´e, Inferring potential landscapes from noisy trajectories of particles within an optical feedback trap, iScience 25 (2022) 104731. doi:10.1016/j.isci.2022.104731.
[13] S. Press´e, I. Sgouralis, Data Modeling for the Sciences, Cambridge University Press, Cambridge, UK, 2023. doi:10.1017/9781009089555.
[14] J. S. Bryan, S. Press´e, Learning continuous potentials from smfret, Biophysical Journal 122 (2023) 433–441. doi:10.1016/j.bpj.2022.11.2947.
[15] N. Masuda, L. E. C. Rocha, A Gillespie algorithm for non-Markovian stochastic processes, SIAM Review 60 (2018) 95–115. doi:10.1137/16M1055876.
[16] B. Munsky, M. Khammash, The finite state projection algorithm for the solution of the chemical master equation, The Journal of Chemical Physics 124 (2006) 044104. doi:10.1063/1.2145882.
[17] P. Pessoa, M. Schweiger, S. Presse, Avoiding matrix exponentials for large transition rate matrices, Journal of Chemical Physics 160 (2024) 094109.
[18] C. Moler, C. Van Loan, Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later, SIAM Review 45 (2003) 3–49. doi:10.1137/S00361445024180.
[19] N. Friedman, L. Cai, X. S. Xie, Linking stochastic dynamics to population distribution: an analytical framework of gene expression., Physical Review Letters 97 (2006) 168302. doi:10.1103/PhysRevLett.97.168302.
[20] V. Shahrezaei, P. S. Swain, Analytical distributions for stochastic gene expression, Proceedings of the National Academy of Sciences 105 (2008) 17256–17261. doi:10.1073/pnas.0803850105.
[21] H. Pendar, T. Platini, R. V. Kulkarni, Exact protein distributions for stochastic models of gene expression using partitioning of poisson processes., Physical Review E 87 (2013) 042720. doi:10.1103/PhysRevE.87.042720.
[22] N. Kumar, T. Platini, R. V. Kulkarni, Exact distributions for stochastic gene expression models with bursting and feedback., Physical Review Letters 113 (2014) 268105. doi:10.1103/PhysRevLett.113.268105.
[23] M. Chen, S. Luo, M. Cao, C. Guo, T. Zhou, J. Zhang, Exact distributions for stochastic gene expression models with arbitrary promoter architecture and translational bursting, Physical Review E 105 (2022) 014405. doi:10.1103/PhysRevE.105.014405.
[24] B. Wu, J. Holehouse, R. Grima, C. Jia, Solving the time-dependent protein distributions for autoregulated bursty gene expression using spectral decomposition, The Journal of Chemical Physics 160 (2024) 074105. doi:10.1063/5.0188455.
[25] C. Jia, Simplification of Markov chains with infinite state space and the mathematical theory of random gene expression bursts., Physical Review E 96 (2017) 032402. doi:10.1103/PhysRevE.96.032402.
[26] F. Jiao, J. Li, T. Liu, Y. Zhu, W. Che, L. Bleris, C. Jia, What can we learn when fitting a simple telegraph model to a complex gene expression model?, PLOS Computational Biology 20 (2024) e1012118. doi:10.1371/journal.pcbi.1012118.
[27] J. Li, S. Wang, W. J. VanDusen, L. D. Schultz, H. A. George, W. K. Herber, H. J. Chae, W. E. Bentley, G. Rao, Green fluorescent protein in Saccharomyces cerevisiae: Real-time studies of the GAL1 promoter, Biotechnology and Bioengineering 70 (2000) 187–196. doi:/10.1002/1097-0290(20001020)70:2¡187::AID-BIT8¿3.0.CO;2-H.
[28] E. Eden, N. Geva-Zatorsky, I. Issaeva, A. Cohen, E. Dekel, T. Danon, L. Cohen, A. Mayo, U. Alon, Proteome half-life dynamics in living human cells, Science 331 (2011) 764–768. doi:10.1126/science.1199784.
[29] E. V. Rusilowicz-Jones, S. Urb´e, M. J. Clague, Protein degradation on the global scale, Molecular Cell 82 (2022) 1414–1423. doi:10.1016/j.molcel.2022.02.027.
[30] J. T. Sauls, D. Li, S. Jun, Adder and a coarse-grained approach to cell size homeostasis in bacteria, Current Opinion in Cell Biology 38 (2016) 38–44. doi:10.1016/j.ceb.2016.02.004.
[31] G. Le Treut, F. Si, D. Li, S. Jun, Quantitative examination of five stochastic cellcycle and cell-size control models for escherichia coli and bacillus subtilis, Frontiers in Microbiology (2021). doi:10.3389/fmicb.2021.721899.
[32] C. Jia, R. Grima, Frequency domain analysis of fluctuations of mRNA and protein copy ers within a cell lineage: Theory and experimental validation, Physical Review X 11 (2021). doi:10.1103/physrevx.11.021032.
[33] A. Torres, S. Cockerell, M. Phillips, G. Bal´azsi, K. Ghosh, Maxcal can infer models from coupled stochastic trajectories of gene expression and cell division, Biophysical Journal 122 (2023) 2623–2635. doi:10.1016/j.bpj.2023.05.017.
[34] Y. Wang, Z. Yu, R. Grima, Z. Cao, Exact solution of a three-stage model of stochastic gene expression including cell-cycle dynamics, The Journal of Chemical Physics 159 (2023). doi:10.1063/5.0173742.
[35] L. Dingeldein, P. Cossio, R. Covino, Simulation-based inference of single-molecule force spectroscopy, Machine Learning Science and Technology 91 (2023) 102988. doi:10.1088/2632-2153/acc8b8.
[36] A. Sukys, K. Ocal, R. Grima, Approximating solutions of the chemical master equation ¨ using neural networks, iScience 25 (2022) 105010. doi:10.1016/j.isci.2022.105010.
[37] L. Dingeldein, P. Cossio, R. Covino, Simulation-based inference of single-molecule experiments, Curr. Opin. Struct. Biol. (2025) 102988doi:10.1016/j.sbi.2025.102988.
[38] K. Cranmer, J. Brehmer, G. Louppe, The frontier of simulation-based inference, Proceedings of the National Academy of Sciences 117 (2020) 30055–30062. doi:10.1073/pnas.1912789117.
[39] M. A. Ramirez Sierra, T. R. Sokolowski, Ai-powered simulation-based inference of a genuinely spatial-stochastic gene regulation model of early mouse embryogenesis, PLOS Computational Biology 20 (2024) e1012473. doi:10.1371/journal.pcbi.1012473.
[40] M. Dax, J. Wildberger, S. Buchholz, S. R. Green, J. H. Macke, B. Sch¨olkopf, Flow matching for scalable simulation-based inference (2023). arXiv:2305.17161, doi:10.48550/arXiv.2305.17161.
[41] D. S. Greenberg, M. Nonnenmacher, J. H. Macke, Automatic posterior transformation for likelihood-free inference, in: International Conference on Machine Learning, 2019. doi:10.48550/arXiv.1905.07488.
[42] G. Papamakarios, I. Murray, Fast ϵ-free inference of simulation models with bayesian conditional density estimation, in: D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, R. Garnett (Eds.), Advances in Neural Information Processing Systems, Vol. 29, Curran Associates, Inc., 2016. doi:10.48550/arXiv.1605.06376.
[43] J. Hermans, A. Delaunoy, F. Rozet, A. Wehenkel, V. Begy, G. Louppe, A crisis in simulation-based inference? beware, your posterior approximations can be unfaithful, Transactions on Machine Learning Research 2022 (2022). doi:/10.48550/arXiv.2110.06581.
[44] M. Falkiewicz, N. Takeishi, I. Shekhzadeh, A. Wehenkel, A. Delaunoy, G. Louppe, A. Kalousis, Calibrating neural simulation-based inference with differentiable coverage probability (2023). arXiv:2310.13402, doi:10.48550/arXiv.2310.13402.
[45] A. Fengler, L. N. Govindarajan, T. Chen, M. J. Frank, Likelihood approximation networks (LANs) for fast inference of simulation models in cognitive neuroscience, eLife 10 (2021) e65074. doi:10.7554/eLife.65074.
[46] T. P. Xiao, Training neural networks using physical equations of motion, Proceedings of the National Academy of Sciences 121 (2024) e241191312. doi:10.1073/pnas.2411913121.
[47] H. Du, C. Krause, V. Mikuni, B. Nachman, I. Pang, D. Shih, Unifying simulation and inference with normalizing flows, Physical Review D 111 (2025) 076004. doi:10.1103/physrevd.111.076004.
[48] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, B. Poole, Score-based generative modeling through stochastic differential equations, in: International Conference on Learning Representations, 2021. doi:10.48550/arXiv.2011.13456.
[49] D. Ghio, Y. Dandi, F. Krzakala, L. Zdeborov´a, Sampling with flows, diffusion, and autoregressive neural networks from a spin-glass perspective, Proceedings of the National Academy of Sciences 121 (2024) e231181012. doi:10.1073/pnas.2311810121.
[50] Q. Jiang, X. Fu, S. Yan, R. Li, W. Du, Z. Cao, F. Qian, R. Grima, Neural network aided approximation and parameter inference of non-Markovian models of gene expression, Nature Communications 12 (2021) 2618. doi:10.1038/s41467-021-22919-1.
[51] C. Liu, J. Wang, Distilling dynamical knowledge from stochastic reaction networks, Proceedings of the National Academy of Sciences 121 (2024) e231742212. doi:10.1073/pnas.2317422121.
[52] E. Nieves, R. Dandekar, C. Rackauckas, Uncertainty quantified discovery of chemical reaction systems via Bayesian scientific machine learning, Frontiers in Systems Biology 4 (2024) 1338518. doi:10.3389/fsysb.2024.1338518.
[53] B. L. Trippe, R. E. Turner, Conditional density estimation with Bayesian normalising flows (2018). arXiv:1802.04908.
[54] G. Papamakarios, D. Sterratt, I. Murray, Sequential neural likelihood: Fast likelihoodfree inference with autoregressive flows, in: K. Chaudhuri, M. Sugiyama (Eds.), Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics, Vol. 89 of Proceedings of Machine Learning Research, 2019, pp. 837–848.
[55] C. Winkler, D. Worrall, E. Hoogeboom, M. Welling, Learning likelihoods with conditional normalizing flows (2023). arXiv:1912.00042.
[56] S. Zhai, R. Zhang, P. Nakkiran, D. Berthelot, J. Gu, H. Zheng, T. Chen, M. A. Bautista, N. Jaitly, J. Susskind, Normalizing flows are capable generative models (2025). URL https://arxiv.org/abs/2412.06329
[57] J. A. Martinez, M. Delvenne, L. Henrion, F. Moreno, S. Telek, C. Dusny, F. Delvigne, Controlling microbial co-culture based on substrate pulsing can lead to stability through differential fitness advantages, PLOS Computational Biology 18 (2022) 1–27. doi:10.1371/journal.pcbi.1010674.
[58] D. Rezende, S. Mohamed, Variational inference with normalizing flows, in: F. Bach, D. Blei (Eds.), Proceedings of the 32nd International Conference on Machine Learning, Vol. 37 of Proceedings of Machine Learning Research, PMLR, Lille, France, 2015, pp. 1530–1538.
[59] G. Papamakarios, T. Pavlakou, I. Murray, Masked autoregressive flow for density estimation, in: I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, R. Garnett (Eds.), Advances in Neural Information Processing Systems, Vol. 30, Curran Associates, Inc., 2017. doi:10.48550/arXiv.1705.07057.
[60] C. Durkan, A. Bekasov, I. Murray, G. Papamakarios, Neural spline flows, in: H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alch´e-Buc, E. Fox, R. Garnett (Eds.), Advances in Neural Information Processing Systems, Vol. 32, Curran Associates, Inc., 2019. doi:10.48550/arXiv.1906.04032.
[61] G. Papamakarios, E. Nalisnick, D. J. Rezende, S. Mohamed, B. Lakshminarayanan, Normalizing flows for probabilistic modeling and inference, Journal of Machine Learning Research. 22 (2021) 57. doi:10.48550/arXiv.1912.02762.
[62] V. Stimper, D. Liu, A. Campbell, V. Berenz, L. Ryll, B. Sch¨olkopf, J. M. Hern´andez-Lobato, normflows: A PyTorch package for normalizing flows, Journal of Open Source Software 8 (2023) 5361. doi:10.21105/joss.05361.
[63] P. Pessoa, M. Schweiger, L. W. Q. Xu, T. Manha, A. Saurabh, J. A. Camarena, S. Press´e, Avoiding subtraction and division of stochastic signals using normalizing flows: Nfdeconvolve (2025). arXiv:2501.08288, doi:10.48550/arXiv.2501.08288.
[64] S. Kullback, R. A. Leibler, On information and sufficiency, The Annals of Mathematical Statistics 22 (1951) 79–86. doi:10.1214/aoms/1177729694.
[65] J. Shore, R. Johnson, Axiomatic derivation of the principle of maximum entropy and the principle of minimum cross-entropy, IEEE Transactions on information theory 26 (1980) 26–37. doi:10.1109/TIT.1980.1056144.
[66] K. Vanslette, Entropic updating of probabilities and density matrices, Entropy 19 (2017) 664. doi:10.3390/e19120664.
[67] P. Pessoa, F. X. Costa, A. Caticha, Entropic dynamics on Gibbs statistical manifolds, Entropy 23 (2021) 494. doi:10.3390/e23050494.
[68] P. Pessoa, Code accompanying "Simulation-trained conditional normalizing flows for likelihood approximation: a case study in stress regulation kinetics in yeast", https://github.com/PessoaP/FC-inference (2025).
[69] F. Jafarpour, C. S. Wright, H. Gudjonson, J. Riebling, E. Dawson, K. Lo, A. Fiebig, S. Crosson, A. R. Dinner, S. Iyer-Biswas, Bridging the timescales of single-cell and population dynamics, Phys. Rev. X 8 (2018) 021007. doi:10.1103/PhysRevX.8.021007.
[70] A. Golubev, Applications and implications of the exponentially modified gamma distribution as a model for time variabilities related to cell proliferation and gene expression, Journal of Theoretical Biology. 393 (2016) 203–217. doi:10.1016/j.jtbi.2015.12.027.
[71] H. X. Chao, R. I. Fakhreddin, H. K. Shimerov, K. M. Kedziora, R. J. Kumar, J. Perez, J. C. Limas, G. D. Grant, J. G. Cook, G. P. Gupta, J. E. Purvis, Evidence that the human cell cycle is a series of uncoupled, memoryless phases, Molecular Systems Biology 15 (2019) e8604. doi:10.15252/msb.20188604.
[72] G. Belluccini, M. L´opez-Garc´ıa, G. Lythe, C. Molina-Par´ıs, Counting generations in birth and death processes with competing Erlang and exponential waiting times, Scientific Reports 12 (2022) 11289. doi:10.1038/s41598-022-14202-0.
[73] H. Haario, E. Saksman, J. Tamminen, An adaptive Metropolis algorithm, Bernoulli 7 (2001) 223 – 242.
[74] S. F. Levy, N. Ziv, M. L. Siegal, Bet hedging in yeast by heterogeneous, age-correlated expression of a stress protectant, PLoS Biol. 10 (2012) e1001325. doi:10.1371/journal.pbio.1001325.
[75] B. M. Zid, E. K. O'Shea, Promoter sequences direct cytoplasmic localization and translation of mRNAs during starvation in yeast, Nature 514 (2014) 117–121. doi:10.1038/nature13578.
[76] L. Henrion, J. A. Martinez, V. Vandenbroucke, M. Delvenne, S. Telek, A. Zicler, A. Gr¨unberger, F. Delvigne, Fitness cost associated with cell phenotypic switching drives population diversification dynamics and controllability, Nat. Commun. 14 (2023) 6128. doi:10.1038/s41467-023-41917-z.
[77] A. I. Konokhova, A. A. Gelash, M. A. Yurkin, A. V. Chernyshev, V. P. Maltsev, Highprecision characterization of individual e. coli cell morphology by scanning flow cytometry, Cytometry Part A 83A (2013) 568–575. doi:10.1002/cyto.a.22294.
[78] BD Biosciences, BD Accuri C6 Plus System User's Guide (2019). URL https://www.bdbiosciences.com/content/dam/bdb/marketing-documents/ BD Accuri C6Plus System Users Guide 2019.pdf
[79] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, E. Duchesnay, scikit-learn: Machine Learning in Python (2011). doi:10.48550/arXiv.1201.0490.
[80] Y. Saeys, S. Van Gassen, B. N. Lambrecht, Computational flow cytometry: helping to make sense of high-dimensional immunology data, Nature Reviews Immunology 16 (2016) 449–462. doi:10.1038/nri.2016.56.
[81] A. K. Kimball, L. M. Oko, B. L. Bullock, R. A. Nemenoff, L. F. van Dyk, E. T. Clambey, A beginner's guide to analyzing and visualizing mass cytometry data, The Journal of Immunology 200 (2018) 3–22. doi:10.4049/jimmunol.1701494.
[82] S. Reich, S. Weissmann, Fokker–planck particle systems for bayesian inference: Computational approaches, SIAM/ASA Journal on Uncertainty Quantification 9 (2021) 446–482. doi:10.1137/19m1303162.
[83] A. Tabandeh, N. Sharma, L. Iannacone, P. Gardoni, Numerical solution of the fokker–planck equation using physics-based mixture models, Computer Methods in Applied Mechanics and Engineering 399 (2022) 115424. doi:10.1016/j.cma.2022.115424.
[84] O. Abdin, P. M. Kim, Direct conformational sampling from peptide energy landscapes through hypernetwork-conditioned diffusion, Nature Machine Intelligence 6 (2024) 775–786. doi:10.1038/s42256-024-00860-4.

Supplemental Information

S1 Why we can assume we know the average cell division time

This Supplementary Information section explains why we assume the mean cell division time T is known in Sec. 3.2 and Sec. 3.3.

In the chemostat experiments for S. cerevisiae, as shown in Sec. 3.4, the average division rate fueled by continuously supplied nutrients matches the dilution rate D of the chemostat. This is achieved through adaptation. Briefly, adaptation occurs by virtue of constraints on growth: if cells divide too slowly, the population is washed out; if they divide too quickly, nutrient depletion or crowding will suppress further growth. This balance can be captured by a simple ordinary differential equation describing the population size N(t):

$$\frac{dN}{dt} = \mu N - DN,\tag{S1}$$

where µ is the specific growth (division) rate of the cells. At steady state, dN dt = 0, we must have µ = D. Since the mean division time, T, is the inverse of the division rate, it follows that T = 1/D.

For this reason, throughout our work, we have adopted units where the mean division time T defines the unit of time, setting T = 1. All results are thus expressed in terms of T. In Sec. 3.4, when analyzing real experimental data, we appropriately convert these dimensionless results back into physical units (hours) by dividing all rates (which have units of inverse time) by D.

S2 Autofluorescence calibration

The distribution of autofluorescence intensity, Iaf, introduced in Sec. 3.3, appears in Eq. (22), but its sampling procedure was not previously specified.

Among the tested conditions, 0.33 h −1 represents the highest dilution rate at which a stable population could be maintained. Since the fluorescent reporter is driven by a stressresponse gene, this condition—with continuous nutrient availability and minimal stress—is assumed to reflect pure autofluorescence. It serves as the reference for autofluorescence calibration, as detailed in Supplemental Information S2.

The distribution for Iaf is modeled using a normalizing flow trained to approximate the distribution of intensities observed at 0.33 h −1 . This normalizing flow uses the same architecture as those presented in the main text but is trained unconditionally—that is, it learns a probability distribution without any conditioning variables. Further implementation details are available on GitHub [68].

The data from the 0.33 h ⁻¹ dilution condition, like the others, was preprocessed as described in SI S5.

Figure S1: Histogram of fluorescence intensities at dilution rate 0.33 h −1 , representing a low-stress control condition dominated by autofluorescence, equivalent to the "control" histogram in Fig. 6a. This distribution was used to train the normalizing flow model defining the autofluorescence intensity distribution Iaf used in (22).

S3 Calibration of fluorescence intensity

To a attest the relationship of the measured florescence in the FC and the number of fluorescent proteins in the cell we used a set of AcGFP1/EGFP flow cytometer calibration beads (Takara Bio, San Jose,CA,cat:632594). This set consists of a mixture of six discrete bead populations having different fluorescent intensities. Each bead population fluorescent intensity is calibrated to a specific molecular equivalent of soluble fluorofor (MESF). Therefore this beads populations were used as the stardards for our intended calibration.

The observed fluorescence intensity I^g is assumed to be proportional to the number of fluorescent proteins G in the cell. As described in the main text, we model this relationship as

$$I_g \sim \text{Normal}(\xi G, \xi G),$$ (S2)

where ξ is a calibration factor that relates molecular counts to fluorescence units. This factor was experimentally determined using calibration beads with known fluorophore content.

Specifically, we used a commercial bead mixture containing six distinct bead populations, each with a known equivalent number of GFP molecules, as provided by the manufacturer. The bead sample was measured under the same flow cytometry settings as the experimental data. We then fit a mixture of six Gaussian distributions to the measured intensity histogram, corresponding to the six known bead populations.

From this fit, we extracted the mean fluorescence intensity of each population. A linear regression was then performed between these mean intensity values and the manufacturerspecified equivalent number of fluorophores. The resulting slope gives the calibration factor ξ, which we estimate to be ξ = 0.0533 fluorescence units per fluorophore.

S4 MCMC strategy

MCMC procedure In this section, we describe the Markov Chain Monte Carlo (MCMC) procedure referenced in Secs. 3.3.2 and 3.4. Although the two applications differ in terms of model and data (with the results in Sec. 3.4 corresponding to the real data setting), both implementations follow a common adaptive MCMC framework. We first present the shared structure of the algorithm, followed by additional details specific to each case.

In both cases, the posterior distribution p(θ|data) is approximated using a normalizing flow, as described in the main text. The flow allows for efficient and differentiable evaluation of the likelihood, which is essential for the computation of the log-posterior during sampling.

Initialization. Prior to sampling, we performed a preliminary search to identify a suitable starting point. We sampled 1,000 parameter vectors θ from the prior distribution and evaluated their unnormalized log-posterior. Each sample was then progressively rescaled toward the current best parameter θ ^∗ using

$$\theta_{\text{prop}} = \theta^* + \frac{1}{\alpha} (\theta^{(i)} - \theta^*), \quad \alpha = 1, \dots, 10.$$ (S3)

If a proposed θprop yielded a higher posterior value, it replaced θ ∗ . This procedure concentrated the search near high-probability regions and produced an initial point close to the posterior mode.

Adaptive phase. We then applied a Metropolis-Hastings MCMC algorithm to generate samples from the posterior distribution p(θ | data). At each step, a candidate parameter θprop ∼ q(· | θ) was drawn from a symmetric Gaussian proposal centered at the current value. The proposed value was accepted with probability

$$A(\theta{\text{prop}}, \theta) = \min\left(1, \frac{p(\theta{\text{prop}}|\text{data})}{p(\theta|\text{data})}\right),$$ (S4)

and otherwise the chain remained at θ.

Inpired by Haario [73], to improve sampling efficiency, the proposal covariance matrix Σ was adapted during the early phase of sampling. Initially, Σ was set to a small isotropic matrix. Every 150 iterations, Σ was updated using the empirical covariance of the most recent 200 samples:

$$\Sigma = \frac{2.4^2}{d} \left( \text{Cov}[\theta^{(t)}] + \epsilon I \right), \tag{S5}$$

where d is the dimensionality of θ, and ϵ = 10−⁸ is a regularization term to ensure numerical stability. Adaptation was performed only if the acceptance rate fell outside the target range of 20–50%. Once the acceptance rate remained within this range for 10 consecutive adaptation steps – meaning that Σ was no longer being updated – the proposal distribution was fixed for the remainder of the chain. This adaptive strategy allowed the proposal to better reflect the geometry of the posterior, improving convergence while maintaining detailed balance.

Sampling stage. After adaptation was complete, we continued sampling using the final proposal distribution. From this stabilized phase of the chain, we collected 10,000 samples to approximate the posterior distribution over parameters. Only these post-adaptation samples are shown in the results presented in Figures 5 and 6.

MCMC initialization in Section 3.4. For the MCMC results presented in Section 3.4, we begin by applying an approximate Bayesian computation (ABC) procedure. In this phase, we simulate data under the model for each proposed parameter value and compute a distance to the observed data based on intensity histograms. A procedure similar to the one described above is then applied to generate 10,000 accepted samples. This initial sampling takes approximately four hours but is not reliable: the likelihood estimates for a fixed parameter value often vary by orders of magnitude, and the method tends to identify only a restricted region of parameter space where the observed fluorescence intensities are approximately matched. Nonetheless, we use the mean and standard deviation of these samples to define a Gaussian prior for the subsequent MCMC sampling described above.

S5 Cleaning real dataset

This section describes the data preprocessing pipeline applied to the flow cytometry measurements discussed in Section 3.4 and visualized in Figure 6a. The goal of this cleaning procedure is to isolate a consistent population of cells suitable for downstream analysis by removing artifacts and outliers that arise during acquisition or due to biological heterogeneity.

The dataset includes three channels of interest: FSC-A (forward scatter area), which relates to cell size; FSC-H (forward scatter height), which helps distinguish single cells from aggregates; and FL1-A, a fluorescence channel used to report gene expression.

As a first step, we removed events with clearly invalid measurements. Specifically, we excluded: Events containing missing NaN values, and non-positive values in any of the key channels which correspond to debris or saturated detectors. We also filtered out events where the ratio between FSC-A and FSC-H deviated by more than 30% from 1. This ratio is expected to be close to unity for round, single cells, and deviations typically suggest doublets or irregular morphologies.

To further refine the population, we applied a density-based outlier detection method aimed at identifying and removing events that lie outside the main cluster of the data. After log-transforming the FSC-A, FSC-H, and FL1-A measurements, we standardized the resulting feature space to remove correlations and ensure comparability between variables. Two empirical filters were also applied to focus the analysis on the densest region of the distribution, excluding events with extreme FSC-H/FSC-A ratios or low FSC-A values.

We then trained a one-class support vector machine (SVM) [79] on a representative subset of this core population. The model was used to define a high-density region in feature space, effectively learning a boundary that separates the main population from peripheral or sparse regions. This approach enables consistent removal of rare or anomalous events across samples, without requiring manual gating.

This cleaning procedure was applied independently to each dilution condition to account for differences in cell density and distribution across experimental setups. The results of the filtering and outlier removal steps are illustrated in Fig. S2, where events removed by the algorithm are shown in red and retained events in blue. Each panel corresponds to a different pairwise projection of the log-transformed features used for classification, with each column coming from cells cultivated at the same dilution rate.

Figure S2: Outlier removal procedure applied to flow cytometry data across three dilution conditions. Each row corresponds to a different projection of the FSC-A, FSC-H, and FL1-A features. Blue points indicate events retained for analysis; red points were classified as outliers and removed. The procedure was applied independently to each condition, these are equivalent to the respective "high stress", "low stress", and "control" in Fig. 6.

S6 Why we assume active/inactive state

The model in Sec. 3.3, later applied in Sec. 3.4, assumes that the cell switches between two states: an active state producing protein at rate β, and an inactive state with zero production.

While more complex models (e.g., multiple or non-parametric states) are conceivable, they must be reflected in the observable fluorescence data. Here we argue that it is reasonable to assume no production in the inactive state. In the case where there is a single producing state, as described in Sec. 3.1, the expected steady-state protein number G in Eq. (15) is given by ⟨G⟩ = 2 βT, where T is the cell division time. In chemostats, where division is set by the dilution rate D (see SI S1), this becomes ⟨G⟩ = β D . If cells were always producing, slower dilution (and thus longer division times) would lead to higher average protein numbers and a shift in the fluorescence peak.

Yet, in Fig. 6a, the peak remains fixed across dilution rates, suggesting a true inactive (effectively zero-production) state compensates for longer division times.

This is supported by simulations presented in Fig. S3, where a one-state always-producing model is used. There, protein synthesis on the order of 104h −1 (two orders of magnitude lower than the levels inferred in Sec. 3.4) leads to fluorescence distributions that change significantly with the dilution rate, even after accounting for autofluorescence. The discrepancy between these simulated peaks and the experimental data further supports the presence of a nonproducing state.

Figure S3: Simulation of a one-state model with constant protein production. The protein number (left) and resulting fluorescence distributions (right) for hypothetic case where we had the same autofluorescence distribution, but there were a contant protein production rate of 10⁴h −1 (about two order of magnitude smaller than the results presented in Fig. 6). As the dilution rate D decreases, the peaks of the distributions shift to higher values due to longer accumulation times. This behavior is inconsistent with the experimental data in Fig. 6a, where peak the low stress has a peak . The result supports the assumption of an inactive state with zero production in the real system.

S7 Pre-calibration of θ distribution for real data using ABC

In Section 3.4, we mention that training the neural network on unrepresentative regions of parameter space can bias inference. To mitigate this, we initialize training using a prior informed by a simpler SBI method, Approximate Bayesian Computation (ABC). This method estimates a rough likelihood by comparing histograms of simulated flow cytometry intensity values at various parameters θ to the observed data shown in Fig. 6. We then perform the MCMC sampling procedure detailed in Sec. S4. We use this MCMC sampled values to construct the prior for Sec. 3.4. The specifics of the ABC approximation are provided in Sec. S7.1, while Sec. S7.2 explains why our normalizing flows-based approach yields a more reliable posterior.

S7.1 ABC-based likelihood

To construct the likelihood approximation using ABC, we discretize the observed log fluorescence intensity values, denoted log I, into A = 200 equally spaced bins. These bins define intervals of the form (La, La+1), where L⁰ = min(log Idata), L^A = max(log Idata), ∆L = La+1 − La.

Given a parameter vector ^θ, we simulate ^S = 2¹² samples ^I^s from the model ^p(I|θ), as described in Sec. 3.4. The approximate likelihood of the log-intensity under these parameters is then defined as:

$$p(\log I \mid \theta) \approx \frac{1}{\Delta L} \cdot \frac{1}{S+1} \sum{a=0}^{A-1} \mathbf{1}{\log I \in (La, L{a+1})} \left[ \frac{1}{A} + \sum{s=1}^{S} \mathbf{1}{\log I_s \in (La, L{a+1})} \right], \tag{S6}$$

where 1x∈(La,La+1) is an indicator function that equals 1 if x falls within the interval (La, La+1), and 0 otherwise. This likelihood estimate is proportional to the number of simulated values log I^s that fall into the same bin as the observed value, with a small pseudocount (1/A) added to each bin, avoiding zero likelihood.

Using this approximation, we proceed with the adaptive MCMC sampling scheme described in Sec. S4. The mean and standard deviation of the resulting ABC posterior samples are used to define the prior distribution over parameters θ.

Specifically, we assume a prior in which the logarithm of each parameter in θ follows a normal distribution. The means and standard deviations of this distribution, shown in Table S1, correspond to those of the ABC posterior. To ensure sufficient coverage of the posterior, standard deviations smaller than 0.5 are rounded up to the nearest half-integer. This prior also defines the distribution from which synthetic data are drawn during neural network training, as used in the results presented in Sec. 3.4.

S7.2 Advantages of the normalizing flow-based posterior

The normalizing flow posterior offers two principal advantages over the ABC-based likelihood approach.

Computational efficiency. Once trained, the normalizing flow provides a closedform, differentiable approximation of the posterior, enabling rapid sampling. For example, generating 10⁴ samples with the ABC-based MCMC scheme (Sec. 3.4) required approximately 4 hours. In contrast, the same time budget with the normalizing flow–based posterior yielded

Parameter	Prior Mean	Prior Std Dev
−1 ln [h ] β	13.2	0.915
High stress (T = 8.33h)
−1 ln λact [T ]	−0.2	0.5
−1 ln λina [T ]	2.5	1.192
ln σ	−1.4	0.5
′ Low stress (T = 4.35h)
′−1 ln [T ] λact	−2.2	0.5
′−1 ln [T ] λina	1.4	1.132
ln σ	−2.8	0.714

Table S1: Prior means and standard deviations for model parameters for the real data analysis in Sec. 3.4. Parameters for activation and inactivation rates (λact and λina) are expressed in units of inverse mean cell division time. Under high stress, T = 1/γ = 8.33h with dilution rate 0.12h −1 . Under low stress, T ′ = 1/γ′ = 4.35h with γ ′ = 0.23h −1 .

10 times more samples. This efficiency is especially beneficial when repeated inference or model refinement is required.

2. Improved MCMC mixing. The normalizing flow posterior also enhances MCMC mixing behavior.

To compare the behavior of the ABC-based and normalizing flows posteriors shown in Fig. 6b, we trace the log-posterior values across MCMC iterations in both cases. In these plots, we subtract the maximum log-posterior value such that the maximum is zero, making relative differences more apparent. This transformation does not affect MCMC dynamics but highlights stability and convergence.

For the ABC-based posterior, we observe that after the initialization described in Sec. S4, the chain does reach regions of higher posterior density. However, even after apparent convergence, the log-posterior continues to oscillate with a magnitude on the order of hundreds. Since the posterior is proportional to the exponential of the log-posterior, these fluctuations imply changes in probability mass by factors as large as e ¹⁰⁰ an implausibly large variation in a supposedly stable region of parameter space. This indicates that the ABC approximation introduces significant noise, stemming from the coarse histogram-based likelihood estimate in (S6).

In contrast, the normalizing flows–based posterior shows smoother convergence and stable sampling after burn-in, with log-posterior values remaining within a narrow band around the mode. This suggests that the flow provides a more accurate and numerically stable surrogate for the likelihood, yielding a better-behaved posterior landscape.

Figure S4: Comparison of log-posterior traces across MCMC iterations for the ABC-based approximation (left) and the normalizing flow–based posterior (right). Each trace has been shifted by subtracting the maximum log-posterior value to emphasize differences in convergence behavior. The ABC-based posterior shows large fluctuations even after burn-in, with log-posterior variations on the order of hundreds corresponding to probability changes by factors of e ¹⁰⁰. This instability arises from the histogram-based likelihood approximation. In contrast, the normalizing flow–based posterior converges more cleanly and remains stable around the maximum, indicating better approximation of the true posterior landscape.

pessoasimulationtrainedconditionalnormalizing2025.md 100 KB تاريخچه خام