Before proceeding to discussing the DGA, we start by briefly reviewing how the traditional Gillespie algorithm simulates discrete stochastic processes. For concreteness, in our exposition, we focus on the chemical system shown in Figure 1 consisting of three species, A, B, and C, whose abundances are described by a state vector ${\displaystyle {\displaystyle \mathbf{x}=(x_1,x_2,x_3)}}$. These chemical species can undergo ${\displaystyle N=3}$ chemical reactions, characterized by rate constants, ${\displaystyle {\displaystyle r_i(\mathbf{x})}}$ where ${\displaystyle {\displaystyle i=1,\dots ,3}}$, and a stoichiometric matrix ${\displaystyle S_{i\alpha }}$ whose ith row encodes how the abundance ${\displaystyle {\displaystyle x_{\alpha }}}$ of species ${\displaystyle {\displaystyle \alpha }}$ changes due to reaction ${\displaystyle i}$. Note that in what follows, we will often supress the dependence of the rates ${\displaystyle {\displaystyle r_i(\mathbf{x})}}$ on ${\displaystyle {\displaystyle \mathbf{x}}}$ and simply write ${\displaystyle r_i}$.

In order to simulate such a system, it is helpful to discretize time into small intervals of size ${\displaystyle {\displaystyle \mathrm{\Delta }t\ll 1}}$. The probability that a reaction ${\displaystyle i}$ with rate ${\displaystyle r_i}$ occurs during such an interval is simply ${\displaystyle {\displaystyle r_i\mathrm{\Delta }t}}$. By construction, we choose ${\displaystyle \mathrm{\Delta }t}$ to be small enough that ${\displaystyle {\displaystyle r_i\mathrm{\Delta }t\ll 1}}$ and that the probability that a reaction occurs in any interval ${\displaystyle {\displaystyle \mathrm{\Delta }t}}$ is extremely small and well described by a Poisson process. This means that naively simulating such a process is extremely inefficient because, in most intervals, no reactions will occur.

The Gillespie algorithm circumnavigates this problem by: (1) exploiting the fact that the reactions are independent so that the rate at which any reaction occurs is also described by an independent Poisson process with rate ${\displaystyle {\displaystyle R={\mathrm{\Sigma }}_i{r}_i}}$ and (2) the waiting time distribution ${\displaystyle {\displaystyle p(\tau )}}$ of a Poisson process with rate ${\displaystyle R}$ is the exponential distribution ${\displaystyle {\displaystyle p(\tau )=R{e}^{-R\tau }}}$. The basic steps of the Gillespie algorithm are illustrated in Figure 1.

The simulation begins with the initialization of time and state variables:

where ${\displaystyle t}$ is the simulation time. One then samples the waiting time distribution ${\displaystyle p(\tau )}$ for a reaction to occur to determine when the next the reaction occurs. This is done by drawing a random number ${\displaystyle u}$ from a uniform distribution over ${\displaystyle [0,1]}$ and updating

Note that this time update is a fully differentiable function of the rates ${\displaystyle r_i}$.

In order to determine which of the reactions ${\displaystyle {\displaystyle i^{\mathrm{\prime }}}}$ occurs after a time ${\displaystyle {\displaystyle \tau }}$, we note that probability that reaction ${\displaystyle i}$ occurs is simply given by ${\displaystyle {\displaystyle q_i=r_i/R}}$. Thus, we can simply draw another random number ${\displaystyle {\displaystyle u^{\mathrm{\prime }}}}$ and choose ${\displaystyle i^{\mathrm{\prime }}}$ such that ${\displaystyle i^{\mathrm{\prime }}}$ equals the smallest integer satisfying

The reaction abundances ${\displaystyle \mathbf{x}}$ are then updated using the stoichiometric matrix

Unlike the time update, both the choice of the reaction ${\displaystyle i^{\mathrm{\prime }}}$ and the abundance updates are not differentiable since the choice of the reaction ${\displaystyle i^{\mathrm{\prime }}}$ is a discontinuous function of the parameters ${\displaystyle r_i}$.

In order to make use of modern deep learning techniques and modern automatic differentiation packages, it is necessary to modify the Gillespie algorithm in such as way as to make the choice of reaction index (Equation 2) and abundance updates (Equation 3) differentiable functions of the kinetic parameters. To do so, we rewrite Equation 2 using a sum of Heaviside step function ${\displaystyle {\displaystyle \mathrm{\Theta }(y)}}$ (recall ${\displaystyle {\displaystyle \mathrm{\Theta }(y)=0}}$ if ${\displaystyle {\displaystyle y<0}}$ and ${\displaystyle \mathrm{\Theta }(y)=1}$ if ${\displaystyle y>0)}$:

This formulation of index selection makes clear the source of non-differentiability. The derivative of the ${\displaystyle i^{\mathrm{\prime }}}$ with respect to ${\displaystyle r_i}$ does not exist at the transition points where the Heaviside function jumps (see Figure 1b).

This suggests a natural modification of the Gillespie algorithm to make it differentiable – replacing the Heaviside function ${\displaystyle \mathrm{\Theta }(y)}$ by a sigmoid function of the form

where we have introduced a ‘hyper-parameter’ ${\displaystyle a}$ that controls the steepness of the sigmoid and plays an analogous role to temperature in a Fermi function in statistical mechanics. A larger value of ${\displaystyle a^{-1}}$ results in a steeper slope for the sigmoid functions, thereby more closely approximating the true Heaviside functions which is recovered in the limit ${\displaystyle a\to 0}$ (see Figure 1b). With this replacement, the index selection equation becomes

Note that in making this approximation, our index is no longer an integer, but instead can take on all real values between 0 and ${\displaystyle N}$. However, by making ${\displaystyle a}$ sufficiently small, Equation 6 still serves as a good approximation to the discrete jumps in Equation 4. In general, ${\displaystyle a}$ is a hyperparameter that is chosen to be as small as possible while still ensuring that the gradient of ${\displaystyle i^{\mathrm{\prime }}}$ with respect to the kinetic parameters ${\displaystyle r_i}$ can be calculated numerically with high accuracy. For a detailed discussion, please see Appendix 1—figure 1 and Appendix 1.

Since the index ${\displaystyle {\displaystyle i^{\mathrm{\prime }}}}$ is no longer an integer but a real number, we must also modify the abundance update in Equation 3 to make it fully differentiable. To do this, we start by rewriting Equation 3 using the Kronecker delta ${\displaystyle {\displaystyle {\delta }_{ij}}}$ (where ${\displaystyle {\displaystyle {\delta }_{ij}=1}}$ if ${\displaystyle {\displaystyle i=j}}$ and ${\displaystyle {\displaystyle {\delta }_{ij}=0}}$ if ${\displaystyle {\displaystyle i\ne j}}$) as

Since ${\displaystyle i^{\mathrm{\prime }}}$ is no longer an integer, we can approximate the Kronecker delta ${\displaystyle {\delta }_{i{i}^{\mathrm{\prime }}}}$ by a Gaussian function, to arrive at the approximate update equation

The hyperparameter ${\displaystyle b}$ is generally chosen to be as small as possible while still ensuring numerical stability of gradients (Appendix 1—figure 1). Note by using an abundance update of the form Equation 8, the species abundances ${\displaystyle \mathbf{x}}$ are now real numbers. This is in stark contrast with the exact Gillespie algorithm where the abundance update (Equation 7) ensures that the ${\displaystyle x_{\alpha }}$ are all integers.

The goal of making Gillespie simulations differentiable is to enable the computation of the gradient of a loss function, ${\displaystyle L(\mathit{\theta })}$, with respect to the kinetics parameters ${\displaystyle \mathit{\theta }}$. A loss function quantifies the difference between simulated and desired values for quantities of interest. For example, when employing the DGA in the context of fitting noisy gene expression models, a natural choice for ${\displaystyle L(\mathit{\theta })}$ is the difference between the simulated and experimentally measured moments of mRNA/protein expression (or alternatively, the Kullback–Leibler divergence between the experimental and simulated mRNA/protein expression distributions if full distributions can be measured). When using the DGA to design gene circuits, the loss function can be any function that characterizes the difference between the simulated and desired values of the input–output relation.

The goal of the optimization using the DGA is to find parameters ${\displaystyle \mathit{\theta }}$ that minimize the loss. The basic workflow of a DGA-based optimization is shown in Figure 2. One starts with an initial guess for the parameters ${\displaystyle {\mathit{\theta }}_0}$. One then uses DGA algorithm to simulate the systems and calculate the gradient of the loss function ${\displaystyle {\mathrm{\nabla }}_{\theta }L(\mathit{\theta })}$. One then updates the parameters, moving in the direction of the gradient using gradient descent or more advanced methods such as ADAM (Kingma and Ba, 2014; Mehta et al., 2019), which uses adaptive estimates of the first and second moments of the gradients to speed up convergence to a local minimum of the loss function.

Figure 2

Download asset Open asset

A summary of the DGA is shown in Figure 1. Unsurprisingly, differentiability comes at a price. The foremost of these is that unlike the Gillespie algorithm, the DGA is no longer exact. The DGA replaces the exact discrete stochastic system by an approximate differentiable stochastic system. This is done by allowing both the reaction index and the species abundances to be continuous numbers. Though in theory, the errors introduced by these approximations can be made arbitrarily small by choosing the hyper-parameters ${\displaystyle a}$ and ${\displaystyle b}$ small enough (see Figure 1), in practice, gradients become numerically unstable when ${\displaystyle a}$ and ${\displaystyle b}$ are sufficiently small (see Appendix 1 and Appendix 1—figure 1).

In what follows, we focus almost exclusively on steady-state properties that probe the ‘bulk’, steady-state properties of the stochastic system of interest. We find the DGA works well in this setting. However, we note that the effect of the approximations introduced by the DGA may be pronounced in more complex settings such as the calculation of rare events, modeling of tail-driven processes, or dealing with non-stationary time series.

In order to better understand the DGA in the context of stochastic gene expression, we benchmarked the DGA on a simple two-state promoter model inspired by experiments in Escherichia coli (Jones et al., 2014). This simple model had several advantages that make it well suited for exploring the performance of DGA. These include the ability to analytically calculate mRNA expression distributions and independent experimental measurements of kinetic parameters.

Gene regulation is tightly regulated at the transcriptional level to ensure that genes are expressed at the right time, place, and in the right amount (Phillips et al., 2012). Transcriptional regulation involves various mechanisms, including the binding of transcription factors to specific DNA sequences, the modification of chromatin structure, and the influence of non-coding RNAs, which collectively control the initiation and rate of transcription (Phillips et al., 2012; Sanchez et al., 2013; Phillips et al., 2019). By orchestrating these regulatory mechanisms, cells can respond to internal signals and external environmental changes, maintaining homeostasis and enabling proper development and function.

Here, we focus on a classic two-state promoter gene regulation (Jones et al., 2014). Two-state promoter systems are commonly studied because they provide a simplified yet powerful model for understanding gene regulation dynamics. These systems, characterized by promoters toggling between active and inactive states, offer insights into how genes are turned on or off in response to various stimuli (see Figure 3a). The two-state gene regulation circuit involves the promoter region, where RNA polymerase (RNAP) binds to initiate transcription and synthesize mRNA molecules at a rate ${\displaystyle r}$. A repressor protein can also bind to the operator site at a rate ${\displaystyle {\displaystyle k_{\text{on}}^{\text{R}}}}$ and unbind at a rate ${\displaystyle {\displaystyle k_{\text{off}}^{\text{R}}}}$. When the repressor is bound to the operator, it prevents RNAP from accessing the promoter, effectively turning off transcription. mRNA is also degraded at a rate ${\displaystyle {\displaystyle \gamma }}$. An appealing feature of this model is that both mean mRNA expression and the Fano factor can be calculated analytically and there exist beautiful quantitative measurements of both these quantities (Figure 3b). For this reason, we use this two-state promoters to benchmark the efficacy of DGA below.

Figure 3

Download asset Open asset

We begin by testing the DGA to do forward simulations on the two-state promoter system described above and comparing the results to simulations performed with the exact Gillespie algorithm (see Appendix 2 for simulation details). Figure 4a compares the probability distribution function (PDF) for the steady-state mRNA levels obtained from the DGA (in red) and the exact Gillespie simulation (in blue). The close overlap of these distributions demonstrates that the DGA can accurately replicate the results of the exact Gillespie simulation. This is also shown by the very close match of the first four moments ${\displaystyle ⟨m^n⟩}$ of the mRNA count between the exact Gillespie and the DGA in Figure 4b, though the DGA systematically overestimates these moments. As observed in Figure 4a, the DGA also fails to accurately capture the tails of the underlying PDF. This discrepancy arises because rare events result from very frequent low-probability reaction events where the sigmoid approximation used in the DGA significantly impacts the reaction selection process and, consequently, the final simulation results.

Figure 4

Download asset Open asset

Next, we compare the accuracy of the DGA in simulating mRNA abundance distributions across a range of simulation times (see Figure 4c). The accuracy is quantified by the ratio of the Jensen–Shannon divergence ${\displaystyle \text{JSD}\left(p_{\text{DGA}}||p_{\text{exact}}^{\text{ss}}\right)}$ between the differentiable Gillespie PDF ${\displaystyle p_{\text{DGA}}}$ and the exact steady-state PDF ${\displaystyle p_{\text{exact}}^{\text{ss}}}$, and the entropy ${\displaystyle \text{H}(p_{\text{exact}}^{\text{ss}})}$ of the exact steady-state PDF. For probability distributions ${\displaystyle P}$ and ${\displaystyle Q}$ over the same discrete space ${\displaystyle \mathcal{X}}$, the JSD and H are defined as:

where ${\displaystyle M=\frac{1}{2}(P+Q)}$ and ${\displaystyle {\displaystyle D_{\text{KL}}}}$ denotes the Kullback–Leibler divergence

The ratio ${\displaystyle \frac{\mathrm{J}\mathrm{S}\mathrm{D}}{\mathrm{H}}}$ normalizes divergence by entropy, enabling meaningful comparison across systems. As expected, the ${\displaystyle \frac{\text{JSD}}{\text{H}}}$ ratio decreases with increasing simulation time, indicating convergence toward the steady-state distribution of the exact Gillespie simulation. By ‘steady-state distribution’, we mean the long-term probability distribution of states that the exact Gillespie algorithm approaches after a simulation time of 10⁴. The saturation of the ${\displaystyle \frac{\text{JSD}}{\text{H}}}$ ratio at approximately 0.003 for long simulation times is due to the finite values of ${\displaystyle {\displaystyle a^{-1}}}$ and ${\displaystyle {\displaystyle b^{-1}}}$. In percentage terms, this ratio represents a 0.3% divergence, meaning that the DGA’s approximation introduces only a 0.3% deviation from the exact distribution, relative to the total uncertainty (entropy) in the exact system.

Finally, the bar plot in Figure 4d shows simulation results for the probability of the promoter being in the ‘OFF’ and ‘ON’ states as predicted by the DGA (in red) and the exact Gillespie simulation (in blue). The differentiable Gillespie overestimates the probability of being in the ‘OFF’ state and underestimates the probability of being in the ‘ON’ state. Nonetheless, given the discrete nature of this system, the DGA does a reasonable job of matching the results of the exact simulations.

As we will see below, despite these errors the DGA is able to accurately capture gradient information and hence works remarkably well at gradient-based optimization of loss functions.

In many applications, one often wants to estimate kinetic parameters from experimental measurements of a stochastic system (Tian et al., 2007; Munsky et al., 2009; Komorowski et al., 2009; Villaverde et al., 2019). For example, in the context of gene expression, biologists are often interested in understanding biophysical parameters such as the rate at which promoters switch between states or a transcription factor unbinds from DNA. However, estimating kinetic parameters in stochastic systems poses numerous challenges because the vast majority of methods for parameter estimation are designed with deterministic systems in mind. Moreover, it is often difficult to analytically calculate likelihood functions making it difficult to perform statistical inference. One attractive method for addressing these difficulties is to combine differentiable Gillespie simulations with gradient-based optimization methods. By choosing kinetic parameters that minimize the difference between simulations and experiments as measured by a loss function, one can quickly and efficiently estimate kinetic parameters and error bars.

To use the DGA for parameter estimation, we start by defining a loss function ${\displaystyle L(\mathit{\theta })}$ that measures the discrepancy between simulations and experiments. In the context of the two-state promoter model (Figure 3), a natural choice of loss function is the square error between the simulated and experimentally measured mean and standard deviations of the steady-state mRNA distributions:

where ${\displaystyle \widehat{⟨m⟩}}$ and ${\displaystyle {\hat{\sigma }}_m}$ denote the mean and standard deviation obtained from DGA simulations, and ${\displaystyle ⟨m⟩}$ and ${\displaystyle {\sigma }_m}$ are the experimentally measured values of the same quantities. Having specified the loss function and parameters, we then use the gradient-based optimization to minimize the loss and find the optimal parameters ${\displaystyle \hat{\mathit{\theta }}}$ (see Figure 2). Note that in general the solution to the optimization problem need not be unique (see below).

Given a set of learned parameters ${\displaystyle \hat{\mathit{\theta }}}$ that minimize ${\displaystyle L(\mathit{\theta })}$, one would also ideally like to assign a confidence interval (CI) to this estimate that reflect how constrained these parameters are. One natural way to achieve this is by examining the curvature of the loss function as the parameter ${\displaystyle {\theta }_i}$ varies around its minimum value, ${\displaystyle {\theta }_i^{\text{min}}}$. Motivated by this, we define the 95% CIs for parameter ${\displaystyle {\theta }_i}$ by:

where

and ${\displaystyle {\displaystyle L({\theta }_i^{\text{min}}-\delta )=L({\theta }_i^{\text{min}}+1.96{\delta }_{{\theta }_i})}}$. A detailed explanation of how to numerically estimate the CIs is given in Appendix 3.

One shortcoming of Equation 13 is that it treats each parameter in isolation and ignores correlations between parameters. On a technical level, this is reflected in the observation that the CIs only know about the diagonal elements of the full Hessian ${\displaystyle {\mathrm{\partial }}_{ij}^{2}L(\mathit{\theta })}$. This shortcoming is especially glaring when there are many sets of parameters that all optimize the loss function (Einav et al., 2018; Razo-Mejia et al., 2018). As discuss below, this is often the case in many stochastic systems including the two-state promoter architecture in Figure 3. For this reason, it is often useful to make two dimensional plots of the loss function ${\displaystyle L(\mathit{\theta })}$. To do so, for each pair of parameters, we simply sample the parameters around their optimal value and forward simulate to calculate the loss function ${\displaystyle L(\mathit{\theta })}$. We then use this simulations to create two-dimensional heat maps of the loss function. This allows us to identify ‘soft directions’ in parameter space, where the loss function ${\displaystyle L(\mathit{\theta })}$ changes slowly, indicating weak sensitivity to specific parameter combinations.

Before proceeding to experiments, we start by benchmarking the DGA’s ability to perform parameter estimation on synthetic data generated using the two-state promoter model shown in Figure 3. This model nominally has four independent kinetic parameters: the rate at which repressors bind the promoter, ${\displaystyle k_{\text{on}}^{\text{R}}}$; the rate at which the repressor unbinds from the promoter, ${\displaystyle k_{\text{off}}^{\text{R}}}$; the rate at which mRNA is produced, ${\displaystyle r}$; and the rate at which mRNA degrades, ${\displaystyle \gamma }$. Since we are only concerned with steady-state properties of the mRNA distribution, we choose to measure time in units of the off rate and set ${\displaystyle k_{\text{off}}^{\text{R}}=1}$ in everything that follows. In Appendix 4, we make use of exact analytical results for ${\displaystyle ⟨m⟩}$ and ${\displaystyle {\sigma }_m}$ to show that the solution to the optimization problem specified by loss function in Equation 11 is degenerate – there are many combinations of the three parameters ${\displaystyle \{k_{\text{on}}^{\text{R}},r,\gamma \}}$ that all optimize ${\displaystyle L(\mathit{\theta })}$. On the other hand, if one fixes the mRNA degradation rate ${\displaystyle \gamma }$, this degeneracy is lifted and there is a unique solution to the optimization problem for the two parameters ${\displaystyle \{k_{\text{on}}^{\text{R}},r\}}$. We discuss both these cases below.

To generate synthetic data, we randomly sample the three parameters: ${\displaystyle k_{\text{on}}^{\text{R}}}$, ${\displaystyle r}$, and ${\displaystyle \gamma }$ within the range ${\displaystyle [0.1,10]}$, while keeping ${\displaystyle k_{\text{off}}^{\text{R}}}$ fixed at 1. In total, we generate 20 different sets of random parameters. We then perform exact Gillespie simulations for each set of parameters. Using these simulations, we obtain the mean ${\displaystyle ⟨m⟩}$ and standard deviation ${\displaystyle {\sigma }_m}$ of the mRNA levels, which are then used as input to the loss function in Equation 11. We then use the DGA to estimate the parameters using the procedure outlined above and compare the resulting predictions with ground truth values for simulations.

We begin by considering the case where the mRNA degradation rate ${\displaystyle \gamma }$ is known and the goal is to estimate the two other parameters: the repressor binding rate ${\displaystyle k_{\text{on}}^{\text{R}}}$ and the mRNA production rate ${\displaystyle r}$. As discussed above, in this case, the loss function in Equation 11 has a unique minima, considerably simplifying the inference task. Figure 5a shows a scatter plot of the learned and the true parameter values for wide variety of choices of ${\displaystyle \gamma }$. As can be seen, there is very good agreement between the true parameters and learned parameters. Figure 5c shows that even when the true and learned parameters differ, the DGA can predict the mean ${\displaystyle ⟨m⟩}$ and standard deviation ${\displaystyle {\sigma }_m}$ of the steady-state mRNA distribution almost perfectly (see Appendix 5 for discussion of how error bars were estimated). To better understand this, we selected a set of learned parameters: ${\displaystyle k_{\text{on}}^{\text{R}}=0.87}$, ${\displaystyle r=3.83}$, and ${\displaystyle \gamma =2.43}$. We then plotted the loss function in the neighborhood of these parameters (Figure 5b). As can be seen, the loss function around the true parameters is quite flat and the learned parameters live at the edge of this flat region. The flatness of the loss function reflects the fact that the mean and standard deviation of the mRNA distribution depend weakly on the kinetic parameters.

Figure 5

Download asset Open asset

We now estimate parameters for the two-state promoter model when all three parameters ${\displaystyle k_{\text{on}}^{\text{R}}}$, ${\displaystyle r}$, and ${\displaystyle \gamma }$ are unknown. As discussed above, in this case, there are many sets of parameters that all minimize the loss function in Equation 11. Figure 6a shows a comparison between the learned and true parameters along with a heat map of the loss function for one set of synthetic parameters (Figure 6b). As can be seen in the plots, though the true parameters and learned parameter values differ significantly, they do so along ‘sloppy’ directions where loss function is flat. Consistent with this intuition, we performed simulations comparing the mean ${\displaystyle \widehat{⟨m⟩}}$ and standard deviation ${\displaystyle {\hat{\sigma }}_m}$ of the steady-state mRNA levels using the true and learned parameters and found near-perfect agreement across all of the synthetic data (Figure 6c).

Figure 6

Download asset Open asset

In the previous section, we demonstrated that our DGA can effectively obtain parameters for synthetic data. However, real experimental data often contains noise and variability, which can complicate the parameter estimation process. To test the DGA in this more difficult setting, we reanalyze experiments by Jones et al., 2014 which measured how mRNA expression changes in a system well described by the two-state gene expression model in Figure 3. In these experiments, two constitutive promoters lacUD5 and 5DL1 (with different transcription rates ${\displaystyle r}$) were placed under the control of a LacI repressor through the insertion of a LacI binding site. By systematically varying LacI concentrations, the authors were able to adjust the repressor binding rate ${\displaystyle k_{\text{on}}^{\text{R}}}$. mRNA fluorescence in situ hybridization was employed to measure mRNA expression, providing data on both mean expression levels ${\displaystyle ⟨m⟩}$ and the variability as quantified by the Fano factor ${\displaystyle f={\sigma }_m^2/⟨m⟩}$ for both promoters (see Figure 3b).

Given a set of measurements of the mean and Fano factor ${\displaystyle \{⟨m{⟩}_i,f_i^m\}}$ for a promoter (lacUD5 and 5DL1), we construct a loss function of the form:

where ${\displaystyle i}$ runs over data points (each with a different lac repressor concentration) and ${\displaystyle \widehat{⟨m{⟩}_i}}$ and ${\displaystyle {\hat{\sigma }}_i^m}$ are the mean and standard deviation obtained from a sample of DGA simulations. This loss function is chosen because, at its minimum, ${\displaystyle \widehat{⟨m{⟩}_i}=⟨m{⟩}_i}$ and ${\displaystyle {\hat{\sigma }}_i^m=\sqrt{f_i^m⟨m{⟩}_i}}$ for all ${\displaystyle i}$. As above, we set ${\displaystyle k_{\text{off}}^{\text{R}}=1}$, and focus on estimating the other three parameters ${\displaystyle \{r,\gamma ,k_{\text{on}}^{\text{R}}\}}$. When performing our gradient-based optimization, we assume that the transcription rate ${\displaystyle r}$ and the mRNA degradation rate ${\displaystyle \gamma }$ are the same for all data points ${\displaystyle i}$, while allowing ${\displaystyle k_{\text{on}}^{\text{R}}}$ to vary across data points ${\displaystyle i}$. This reflects the fact that ${\displaystyle k_{\text{on}}^{\text{R}}}$ is a function of the lac repressor concentration which, by design, is varied across data points (see Appendix 6 for details on how this optimization is implemented and calculation of error bars).

The results of this procedure are summarized in Figure 7. We find that for the lacUD5 promoter ${\displaystyle \hat{r}=90.25}$, ${\displaystyle \hat{\gamma }=6.20}$ and that ${\displaystyle {\hat{k}}_{\text{on}}^{\text{R}}}$ varies from a minimum value of 0.18 to a maximum value of 99.0. For the 5DL1 promoters ${\displaystyle \hat{r}=87.48}$ and ${\displaystyle \hat{\gamma }=9.80}$ and ${\displaystyle {\hat{k}}_{\text{on}}^{\text{R}}}$ varies between 3.64 and 99.0. Recall that we have normalized all rates to the repressor unbinding rate ${\displaystyle k_{\text{off}}^{\text{R}}=1}$. These values indicate that mRNA transcription occurs much faster compared to the unbinding of the repressor, suggesting that once the promoter is in an active state, it produces mRNA rapidly. The relatively high mRNA degradation rates indicate a mechanism for fine-tuning gene expression levels, ensuring that mRNA does not persist too long in the cell, which could otherwise lead to prolonged expression even after promoter deactivation.

Figure 7

Download asset Open asset

As expected, the repressor binding rates decrease with the mean mRNA level (see Figure 7c). The broad range of repressor binding rates shows that the system can adjust its sensitivity to repressor concentration, allowing for both tight repression and rapid activation depending on the cellular context.

Figure 7a shows a comparison between the predictions of the DGA (solid curves) and the experimental data (squares) for mean mRNA levels and the Fano factor . The theoretical curves are obtained by using analytical expression for and from Gillespie, 2007 with parameters estimated from the DGA. We find that for the lacUD5 and the 5DL1 promoters, the mean percentage errors for predictions of the Fano factor are 25% and 28%, respectively (see Appendix 6).

An appealing feature of Jones et al., 2014 is that the authors performed independent experiments to directly measure the normalized transcription rate ${\displaystyle {\displaystyle r/\gamma }}$ (namely the ratio of the transcription rate and the mRNA degradation rate). This allows us to compare the DGA predictions for these parameters to ground truth measurements of kinetic parameters. In Figure 7b, the predictions of the DGA agree remarkably well for both the lacUD5 and 5DL1 promoters.

Another interesting application of the DGA is to design stochastic chemical or biological networks that exhibit a particular behavior. In many cases, this design problem can be reformulated as identifying choices of parameter that give rise to a desired behavior. Here, we show that the DGA is ideally suited for such a task. We focus on designing the input–output relation of a four state promoter model of gene regulation (Lammers et al., 2023). We have chosen this more complex promoter architecture because, unlike the two-state promoter model analyzed above, it allows for nonequilibrium currents. In making this choice, we are inspired by numerous recent works have investigated how cells can tune kinetic parameters to operate out of equilibrium in order to achieve increased sharpness/sensitivity (Nicholson and Gingrich, 2023Lammers et al., 2023; Zoller et al., 2022; Wong and Gunawardena, 2020; Dixit et al., 2024).

We focus on designing the steady-state input–output relationship of the four-state promoter model of gene regulation model shown in Figure 8a; Lammers et al., 2023. The locus can be in either an ‘ON’ state where mRNA is transcribed at a rate ${\displaystyle r}$ or an ‘OFF’ state where the locus is closed and there is no transcription. In addition, a transcription factor (assumed to be an activator) with concentration ${\displaystyle [c]}$ can bind to the locus with a concentration dependent rate ${\displaystyle [c]k_b}$ in the ‘OFF’ state and a rate ${\displaystyle [c]{\eta }_{ba}{k}_b}$ in the ‘ON’ rate. The activator can also unbind at a rate ${\displaystyle k_u}$ in the ‘OFF’ state and a rate ${\displaystyle {\eta }_{ua}{k}_u}$ in the ‘ON’ state. The average mRNA production rate (averaged over many samples) in this model is given by

where ${\displaystyle {\displaystyle {\pi }_s}}$ (${\displaystyle {\displaystyle s=2,3}}$) is the steady-state probability of finding the system in each of the ‘ON’ states (see Figure 8a).

Figure 8

Download asset Open asset

Such promoter architectures are often studied in the context of protein gradient-based development (Lammers et al., 2023; Estrada et al., 2016; Owen and Horowitz, 2023). One well-known example of such a gradient is the dorsal protein gradient in Drosophila, which plays a crucial role in determining the spatial boundaries of gene expression domains during early embryonic development. In this context, the sharpness of the response as a function of activator concentration is a critical aspect. High sharpness ensures that the transition between different gene expression domains occurs over a very narrow region, leading to well-defined and precise boundaries. Inspired by this, our objective is to determine the parameters such that the variation in ${\displaystyle ⟨\overline{r}⟩}$ as a function of the activator concentration ${\displaystyle [c]}$ follows a desired response. We consider the two target responses (shown in Figure 8b) of differing sharpness, which following Lammers et al., 2023 we quantify as ${\displaystyle max\left(\frac{\mathrm{\partial }⟨\overline{r}⟩}{\mathrm{\partial }[c]}[c]\right)}$. For simplicity, we use sixth-degree polynomials to model the input–output functions, with the x-axis plotted on a logarithmic scale. We note that our results do not depend on this choice and any other functional form works equally well.

In order to use the DGA to learn a desired input–output relation, we must specify a loss function that quantifies the discrepancy between the desired and actual responses of the promoter network. To construct such a loss function, we begin by discretizing the activator concentration into ${\displaystyle N=10}$ logarithmically spaced points, ${\displaystyle [c{]}_i}$, where ${\displaystyle i=1,2,\dots ,N}$. For each ${\displaystyle [c{]}_i}$, we denote the corresponding average mRNA production rate ${\displaystyle ⟨\overline{r}{⟩}_i}$ (see Equation 15). After discretization, the loss function is simply the square error between the desired response, ${\displaystyle ⟨\overline{r}{⟩}_i}$, and the current response, ${\displaystyle ⟨\hat{\overline{r}}(\mathit{\theta }){⟩}_i}$, of the circuit

where ${\displaystyle ⟨\hat{\overline{r}}(\mathit{\theta }){⟩}_i}$ denotes the predicted average mRNA production rates obtained from the DGA simulations given the current parameters ${\displaystyle \mathit{\theta }}$. To compute ${\displaystyle ⟨\hat{\overline{r}}{⟩}_i}$ for a concentration ${\displaystyle [c{]}_i}$, we perform ${\displaystyle n=600}$ DGA simulations (indexed by capital letters ${\displaystyle A=1,\dots ,n}$) using the DGA and use these simulations to calculate the fraction of time spent in transcriptionally active states (states ${\displaystyle s=2}$ and ${\displaystyle s=3}$ in Figure 8a). If we denote the fraction of time spent in state ${\displaystyle s}$ in simulation ${\displaystyle A}$ by ${\displaystyle w_s^A}$, then we can calculate the probability ${\displaystyle {\pi }_s}$ of being in state ${\displaystyle s}$ by

and use Equation 15 to calculate ${\displaystyle ⟨\hat{\overline{r}}(\mathit{\theta }){⟩}_i}$

As before, we optimize this loss using gradient descent (see Figure 2). We assume that the transcription rate ${\displaystyle r}$ is known (this just corresponds to an overall scaling of mRNA numbers). Since we are concerned only with steady-state properties, we fix the activator binding rate to a constant value, ${\displaystyle k_b=0.02}$. This is equivalent to measuring time in units of ${\displaystyle k_b^{-1}}$. We then use gradient descent to optimize the remaining seven parameters governing transitions between promoter states.

Figure 8b shows a comparison between the desired and learned input–output relations. This is good agreement between the learned and desired responses, showing that the DGA is able to design dose–response curves with different sensitivities and maximal values. Figure 8c shows the learned parameters for both response curves. Notably, the degree of activation resulting from transcription factor binding, denoted by ${\displaystyle {\eta }_{ab}}$, is substantially higher for the sharper response (Response-2). In contrast, the influence on transcription factor binding due to activation, represented by ${\displaystyle {\eta }_{ba}}$, is reduced for the sharper response curve. Additionally, the unbinding rate ${\displaystyle k_u}$ is observed to be lower for the sharper response. However, it is essential to approach these findings with caution, as the parameters are highly interdependent. These interdependencies can be visualized by plotting the loss function around the optimized parameter values. Figure 8e shows two dimensional heat maps of the loss function for Response-2. There are seven free parameters, resulting in a total of 21 possible 2D slices of the loss function within the seven-dimensional loss landscape.

The most striking feature of these plots is the central role played by the parameters ${\displaystyle {\eta }_{ab}}$ and ${\displaystyle {\eta }_{ua}}$ which must both be high, suggesting that the sharpness in Response-2 may result from creating a high-flux nonequilibrium cycle through the four promoter states (see Figure 8a). This observation is consistent with recent works suggesting that creating such nonequilibrium kinetics represents a general design principle for engineering sharp responses (Lammers et al., 2023; Zoller et al., 2022; Wong and Gunawardena, 2020; Dixit et al., 2024). To better understand if this is indeed what is happening, we quantified the energy dissipation per unit time (power consumption), ${\displaystyle \mathrm{\Phi }}$, in the nonequilibrium circuit. The energetic cost of operating biochemical networks can be quantified using ideas from nonequilibrium thermodynamics using a generalized Ohm’s law of the form (Lammers et al., 2023; Qian, 2007; Mehta and Schwab, 2012; Lan et al., 2012; Lang et al., 2014; Mehta et al., 2016)

where we have defined a nonequilibrium drive

and the nonequilibrium flux

where ${\displaystyle {\pi }_0}$ and ${\displaystyle {\pi }_1}$ are the probabilities of finding the system in state 0 and 1, respectively. Figure 8d shows a comparison between energy consumption and sharpness of the two learned circuits. Consistent with the results of Lammers et al., 2023, we find that the sharper response curve is achieved by consuming more energy.

In this section, we explore the tradeoffs involved in tuning the hyperparameters ${\displaystyle a}$ and ${\displaystyle b}$ in the DGA. These hyperparameters are crucial for balancing the accuracy and numerical stability of the DGA in approximating the exact Gillespie algorithm.

The hyperparameter ${\displaystyle a^{-1}}$ controls the steepness of the sigmoid function used to approximate the Heaviside step function in reaction selection. Similarly, ${\displaystyle b^{-1}}$ determines the sharpness of the Gaussian function used to approximate the Kronecker delta function in abundance updates. A larger value of ${\displaystyle a^{-1}}$ or ${\displaystyle b^{-1}}$ results in a steeper sigmoid or Gaussian function, thus more closely approximating the discrete functions in the exact Gillespie algorithm.

To assess the impact of these hyperparameters on the accuracy of the DGA, we measure the ratio of the Jensen–Shannon divergence between the DGA-generated PDF ${\displaystyle p_{\text{DGA}}}$ and the exact PDF ${\displaystyle p_{\text{exact}}}$, normalized by the entropy ${\displaystyle \text{H}(p_{\text{exact}})}$ of the exact steady-state PDF (see Equation 9). This ratio, ${\displaystyle \frac{\text{JSD}(p_{\text{DGA}}\Vert p_{\text{exact}})}{\text{H}(p_{\text{exact}})}}$, provides a measure of how closely the DGA approximates the exact Gillespie algorithm.

Appendix 1—figure 1 shows the ratio ${\displaystyle \frac{\text{JSD}}{\text{H}}}$ as a function of the sharpness parameters ${\displaystyle a^{-1}}$ and ${\displaystyle b^{-1}}$. In panel (a), ${\displaystyle b^{-1}}$ is fixed at 20, and ${\displaystyle a^{-1}}$ is varied. In panel (b), ${\displaystyle a^{-1}}$ is fixed at 200, and ${\displaystyle b^{-1}}$ is varied. Some key insights can be drawn from these plots.

First, as ${\displaystyle a^{-1}}$ or ${\displaystyle b^{-1}}$ increases, the ratio ${\displaystyle \frac{\text{JSD}}{\text{H}}}$ decreases, indicating that the DGA’s approximation becomes more accurate. This is because steeper sigmoid and Gaussian functions better mimics the discrete steps of the exact Gillespie algorithm. Interestingly, while the ratio decreases for both parameters, ${\displaystyle a^{-1}}$ plateaus at high values, whereas ${\displaystyle b^{-1}}$ rises at high values. This plateau occurs because the sigmoid function used for reaction selection becomes so steep that it effectively becomes a step function, beyond which further steepening has negligible impact. As ${\displaystyle b^{-1}}$ becomes very large, the width of the Gaussian function used for abundance updates becomes extremely narrow. In such a scenario, it becomes increasingly improbable for the chosen reaction index to fall within this narrow width, especially because the reaction indices are not exact integers but are instead near-integer continuous values. Therefore, as ${\displaystyle b^{-1}}$ becomes very large, the discrepancy between the DGA-generated probabilities and the exact probabilities widens, causing the ratio ${\displaystyle \frac{\text{JSD}}{\text{H}}}$ to increase.

Appendix 1—figure 1

Download asset Open asset

Numerical stability for gradient computation is crucial. Therefore, it is important to examine how the gradient behaves as a function of the hyperparameters. Panels (c) and (d) of Appendix 1—figure 1 provide insights into this behavior. The plots show the gradient ${\displaystyle \mathrm{\nabla }{L}_r}$ of the loss function ${\displaystyle L}$ with respect to the parameter ${\displaystyle r}$ near the true parameter values. The gradients are computed for the loss function in Equation 11 with ${\displaystyle ⟨m⟩}$ and ${\displaystyle {\sigma }_m}$ equal to 8 and 2.5, respectively, for the parameter values ${\displaystyle k_{\text{on}}^{\text{R}}=0.5}$, ${\displaystyle k_{\text{off}}^{\text{R}}=1}$, ${\displaystyle r=10}$, and ${\displaystyle \gamma =1}$. With these parameter values, the true mean and standard deviation are equal to 6.67 and 3.94, respectively.

As ${\displaystyle a^{-1}}$ or ${\displaystyle b^{-1}}$ increases, the gradients become more accurate, but their numerical stability can be compromised. This is evidenced by the increased variability and erratic behavior in the gradients at very high sharpness values (see Appendix 1—figure 1). Hence, very large values of ${\displaystyle a^{-1}}$ or ${\displaystyle b^{-1}}$ lead to oscillations and convergence issues, highlighting the need for a balance between accuracy and stability.

The tradeoff between accuracy and numerical stability is evident. This necessitates careful tuning of ${\displaystyle a}$ and ${\displaystyle b}$ to ensure stable and efficient optimization. In practice, this involves selecting values that provide sufficient approximation quality without compromising the stability of the gradients.

This section explains the implementation of the DGA used to simulate the two-state promoter model. The implementation can be adapted to any model where the stoichiometric matrix and propensities are known. The model involves promoter state switching and mRNA production/degradation (Figure 3a), and the algorithm is designed to handle these sub-processes effectively.

The stoichiometric matrix is a key component in modeling state changes due to reactions. Each row represents a reaction, and each column corresponds to a state variable (promoter state and mRNA level). The matrix for the two-state promoter model includes:

• Reaction 1: Promoter state transitions from the OFF state (–1) to the ON state (+1).

• Reaction 2: Production of mRNA.

• Reaction 3: Promoter state transitions from the ON state (+1) to the OFF state (–1).

• Reaction 4: Degradation of mRNA.

The stoichiometric matrix ${\displaystyle S}$ for this model is:

In this matrix, the rows correspond to the reactions listed above, and the columns represent the promoter state and mRNA level, respectively.

The levels vector is a two-dimensional vector where the first element represents the promoter state (–1 or +1) and the second element represents the mRNA number ${\displaystyle m}$. Propensities are the rates at which reactions occur, given the current state of the system. In our PyTorch implementation, the propensities are calculated using the following expressions:

propensities = torch.stack([
   kon * torch.sigmoid(-c * levels[0]),
   # Promoter state switching from –1 to +1
   r * torch.sigmoid(-c * levels[0]),
   # mRNA production
   torch.sigmoid(c * levels[0]),
   # Promoter state switching from +1 to –1
   g * levels [1]
   # mRNA degradation
])

Each propensity corresponds to a different reaction:

• Promoter state switching from –1 to +1: The value of kon * torch.sigmoid(-c * levels[0]) is around ${\displaystyle k_{\text{on}}^{\text{R}}}$ when the levels[0] (promoter state) is around –1 and close to zero when levels[0] is around +1. The sigmoid function ensures a smooth transition, allowing differentiability and preventing abrupt changes in rates. The constant ${\displaystyle c}$ controls the sharpness of the sigmoid function.

• mRNA production: The rate is proportional to ${\displaystyle r}$ and modulated by the same sigmoid function, torch.sigmoid(-c * levels[0]), such that the rate is equal to ${\displaystyle r}$ only when the promoter is in –1 state.

• Promoter state switching from +1 to –1: The rate is set to 1 and modulated by the sigmoid function, torch.sigmoid(c * levels[0]), such that the rate is equal to 1 only when the promoter is in +1 state.

• mRNA degradation: The rate is proportional to the current mRNA level, g * levels [1], reflecting the natural decay of mRNA with rate ${\displaystyle m\gamma }$.

Using the sigmoid function in the propensity calculations is crucial for ensuring smooth and differentiable transitions between states. This smoothness is essential for gradient-based optimization methods, which rely on continuous and differentiable functions to compute gradients effectively. Without the sigmoid function, the propensity rates could change abruptly, leading to numerical instability and difficulties in optimizing the model parameters.

Reaction selection function: We define a function reaction_selection that selects the next reaction to occur based on the transition points and a random number between [0,1]. The function basically implements using Equation 6. The transition points are first calculated from the cumulative sum of reaction rates normalized to the total rate.

Gillespie simulation function: The main gillespie_simulation function uses the previously described functions, each with specific roles, to perform the actual simulation step-by-step, as shown in Figure 1. This function iterates through the number of simulations, updating the system’s state and the propensities of each reaction at each step.

State jump function: We define another function state_jump that calculates the state change vector when a reaction occurs. It uses a Gaussian function to smoothly transition between states based on the selected reaction index and the stoichiometry matrix (see Equation 8).

In this section, we describe the methodology used to estimate the confidence intervals for the parameters using polynomial fitting and numerical techniques. This approach uses the results of multiple simulations to determine the range within which the parameter ${\displaystyle {\theta }_i}$ is likely to lie, based on the curvature of the loss function around its minimum value. Specifically, we perform the following steps:

1. Parameter initialization: Set up and initialize the necessary parameters. This includes defining the number of evaluation points, the number of simulations, the simulation time, and the hyperparameters ${\displaystyle a^{-1}}$, ${\displaystyle b^{-1}}$, and ${\displaystyle c}$.

2. Range generation: For each set of learned parameters ${\displaystyle \hat{\mathit{\theta }}}$, generate a range of values for the parameter of interest ${\displaystyle {\theta }_i}$ while keeping the other parameters fixed at their learned values. Let the learned value of the parameter ${\displaystyle {\theta }_i}$ be ${\displaystyle {\hat{\theta }}_i}$. Then the range of evaluation is approximately ${\displaystyle [0.2{\hat{\theta }}_i,2{\hat{\theta }}_i]}$, depending on the flatness of the loss landscape around its minimum.

3. Simulation and loss calculation: For each value in this range, perform forward DGA simulations to calculate the mean and standard deviation of the results, using which the loss function is computed and stored.

4. Polynomial fitting: Fit a polynomial of degree 6 to the computed loss values across the range of ${\displaystyle {\theta }_i}$. Compute the first and second derivatives of the fitted polynomial to identify the minimum and evaluate the curvature.

5. Identifying minimum: Identify the valid minimum ${\displaystyle {\theta }_i^{\text{min}}}$ of the loss function by solving for the roots of the first derivative and filtering out the points where the second derivative is positive.

6. Curvature and standard deviation calculation: Calculate the curvature of the loss landscape at its minimum using its second derivative at the minimum. An estimation of the standard deviation ${\displaystyle {\delta }_{{\theta }_i}}$ of the parameter ${\displaystyle {\theta }_i}$ is given by:

7. Confidence interval calculation: The loss function is typically asymmetric around its minimum, with the left side usually steeper than the right (see Appendix 3—figure 1). To determine the right error bar, we use ${\displaystyle {\displaystyle {\theta }_i^{\text{min}}+1.96{\delta }_{{\theta }_i}}}$. For the left error bar, we find the point where ${\displaystyle {\displaystyle L({\theta }_i^{\text{min}}-\delta )=L({\theta }_i^{\text{min}}+1.96{\delta }_{{\theta }_i})}}$ (see Appendix 3—figure 1). Therefore, the balanced 95% CI for the parameter ${\displaystyle {\displaystyle {\theta }_i}}$ is given by:

This methodology allows for a robust estimation of the confidence intervals, providing insights into the reliability and precision of the parameter estimates.

Appendix 3—figure 1

Download asset Open asset

We will now demonstrate the existence of degeneracy in the two-state promoter architecture. Setting ${\displaystyle k_{\text{off}}^{\text{R}}=1}$ in the analytical expressions for the mean ${\displaystyle ⟨m⟩}$ and the Fano factor ${\displaystyle f}$ from Jones et al., 2014, we have:

We want to solve Equation D1 for ${\displaystyle r}$ and ${\displaystyle \gamma }$. By isolating ${\displaystyle r}$ in the expression for ${\displaystyle ⟨m⟩}$, we obtain:

Next, we substitute the expression for ${\displaystyle r}$ from Equation D2 into the expression for the ${\displaystyle f}$ in Equation D1:

Expanding and isolating terms involving ${\displaystyle \gamma }$:

Rearranging to solve for ${\displaystyle \gamma }$:

Thus, we obtain:

We substitute back Equation 28 in Equation 23 to obtain ${\displaystyle r}$:

Equations D3 and D4 indicate that for any given ${\displaystyle k_{\text{on}}^{\text{R}}}$, there exists a corresponding pair of values for ${\displaystyle r}$ and ${\displaystyle \gamma }$ that satisfy the equations. Therefore, the solutions are not unique, demonstrating a high degree of degeneracy in the parameter space. This degeneracy arises because multiple combinations of ${\displaystyle {\displaystyle \{r,\gamma ,k_{\text{on}}^{\text{R}}\}}}$ can produce the same observable ${\displaystyle {\displaystyle ⟨m⟩}}$ and ${\displaystyle f}$. The lack of unique solutions is a common issue in parameter estimation for complex systems, where different parameter sets can lead to similar system behaviors.

To resolve the degeneracy, we can fix the value of ${\displaystyle \gamma }$. Now, if the rate ${\displaystyle \gamma }$ is known, we can solve Equation D1 for the other parameters. Starting by rearranging the mean expression from Equation D1:

Substituting Equation D5 in the expression of ${\displaystyle f}$ in Equation D1, we have

Substituting Equation D6 into Equation D5, we obtain:

Equations D6 and D7 provide unique solutions for ${\displaystyle r}$ and ${\displaystyle k_{\text{on}}^{\text{R}}}$ given a known value of ${\displaystyle \gamma }$. By fixing ${\displaystyle \gamma }$, the degeneracy is resolved, and the remaining parameters can be accurately determined.

To assess the reliability of the statistics predicted through DGA-based optimization, we calculate error bars using the following procedure. For the non-degenerate situation, we use the learned parameter values ${\displaystyle {\hat{k}}_{\text{on}}^{\text{R}}}$ and ${\displaystyle \hat{r}}$, and perform forward DGA simulations with many different random seeds. This generates multiple samples of the mean mRNA level ${\displaystyle ⟨m⟩}$ and the standard deviation ${\displaystyle {\sigma }_m}$. These samples provide us with the variability due to the stochastic nature of the simulations. The 95% CIs are then determined using their standard deviations, denoted as ${\displaystyle {\delta }_{⟨m⟩}}$ and ${\displaystyle {\delta }_{{\sigma }_m}}$. For the mean mRNA level, the CI is calculated as:

For the standard deviation of mRNA levels, the CI is calculated as:

To quantitatively assess the goodness of the fit, we define the mean percentage error (MPE) as follows:

where ${\displaystyle {\displaystyle {\hat{f}}_i^m}}$ is the predicted Fano factor for the ith data point. This metric provides a measure of the average discrepancy between the predicted and experimental Fano factors, expressed as a percentage of the experimental values.

The error bars in the ratio ${\displaystyle \hat{r}/\hat{\gamma }}$ are obtained by applying error propagation to the standard deviations ${\displaystyle {\delta }_r}$ and ${\displaystyle {\delta }_{\gamma }}$ of the individual values ${\displaystyle \hat{r}}$ and ${\displaystyle \hat{\gamma }}$. The propagated error is given by:

where ${\displaystyle {\delta }_r}$ and ${\displaystyle {\delta }_{\gamma }}$ are the standard deviations of ${\displaystyle \hat{r}}$ and ${\displaystyle \hat{\gamma }}$, respectively. These standard deviations are obtained using the curvature of the loss function, as discussed earlier.

The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.

This work was supported by NIH NIGMS R35GM119461 to PM and Chan-Zuckerburg Institute Investigator grant to PM. The authors also acknowledge support from the Shared Computing Cluster (SCC) administered by Boston University Research Computing Services. We would also like to thank the Mehta and Kondev groups for useful discussions.

1. Preprint posted: September 25, 2024
2. Sent for peer review: October 10, 2024
3. Reviewed Preprint version 1: January 7, 2025
4. Reviewed Preprint version 2: February 26, 2025
5. Version of Record published: March 17, 2025

You can cite all versions using the DOI https://doi.org/10.7554/eLife.103877. This DOI represents all versions, and will always resolve to the latest one.

© 2025, Rijal and Mehta

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.