Quantifying biochemical reaction rates from static population variability within incompletely observed complex networks

Quantifying biochemical reaction rates within complex cellular processes remains a key challenge of systems biology even as high-throughput single-cell data have become available to characterize snapshots of population variability. That is because complex systems with stochastic and non-linear interactions are difficult to analyze when not all components can be observed simultaneously and systems cannot be followed over time. Instead of using descriptive statistical models, we show that incompletely specified mechanistic models can be used to translate qualitative knowledge of interactions into reaction rate functions from covariability data between pairs of components. This promises to turn a globally intractable problem into a sequence of solvable inference problems to quantify complex interaction networks from incomplete snapshots of their stochastic fluctuations.


Simulation details
All simulations were performed using the standard Doob-Gillespie (1, 2) algorithm until each reaction occurred at least 10 8 times. The resulting high-fidelity distributions from these simulations were considered as the underlying 'true' distribution for each respective system. To obtain the distributions used to test our inference algorithm, we randomly sampled N individual data points (N = 100, 000 if not otherwise specified) from the 'true' distributions corresponding to N independent snapshots of a given system. For this sampling we used the numpy.random.choice function in Python. Fig. 2 of main text. Below we specify the parameters used to simulate the numerical proof-of-principle examples introduced in Figure 2 of the main text. These parameters do not necessarily correspond to realistic cellular networks but were chosen to achieve a distinct set of system dynamics to make the mathematical point that we can successfully infer f (x 2 ) for diverse system dynamics across comparable range for X 2 -variability. All variables are subject to the following reactions

Reaction dynamics for numerical proof-of-principle examples in
with a linear degradation rate dependent on the lifetime τ i of a specific molecule X i and a production rate f i (x) dependent on the state vector x.
As mentioned in the main text the production rate of X 3 was kept the same for all systems: with a lifetime of τ 3 = 1 and the Hill function parameters λ 3 = 80, n 3 = 2, K 3 = 40. The other rate functions were chosen as follows: • Bistable system • Oscillating system Lifetimes: • Noise controlling system ; λ 2 = 3000, n 2 = 10, K 2 = 10 Lifetimes: τ 1 = 50, τ 2 = τ 3 = 1 • Noise enhancing system Lifetimes: Reaction dynamics for example systems in Fig. 3-5 of main text. To illustrate how the inference method depends on time-scales and sampling (Fig. 3 of the main text) we simulated the "noise enhancing system" as defined above while logarithmically varying the life-time τ 3 = 1.5 k with k ranging from −7 to 10. In order to keep system averages x 3 comparable while varying the time-scales we simultaneously adjusted the value of λ 3 such that λ 3 τ 3 = 80 remained constant.
To illustrate that the nature of the degradation rate does not significantly affect the inference (Fig. 5 of the main text) we simulated the above "noise enhancing system" but replaced the degradation reaction of X 3 with with a fixed value of γ = 2.

Modeling measurement errors -Fig. 6 of the main text
In the main text we analyzed how measurement noise affects our inference method by accounting for absolute and relative error terms as well as undercounting of molecules. The following error models (x 2 , x 3 ) → (x 2 , x 3 ) were applied to each sampled data point: • abolute error (x 2 , x 3 ) = (N (x 2 , σ abs ), N (x 3 , σ abs )), • relative error • binominal undercounting where N (µ, σ) refers to a normal distribution with a mean of µ and standard variation of σ. B(n, p) refers to a binomial distribution where p corresponds to the probability that a specific molecule gets detected. All resulting (x 2 , x 3 ) pairs were rounded to integers and pairs with negative molecule numbers were ignored.
When applying our inference method to simulated data with binomial undercounting the resulting rate function depends on the undercounted variable x 2 . The correct rate function was obtained by rescaling the inferred rate function argument with the detection probability p such that f (x 2 ) = f (x 2 /p).

Derivation of the general probability balance relation
We consider a general stochastic system with state vector x = (x 1 , x 2 , x 3 , . . . , x ), where the k th reaction changes levels of component X j by d kj , with reaction propensities r k (x) that are arbitrarily non-linear functions of the state vector. The chemical master equation is then given by Considering one component of interest X i and summing over all other components X j = X i gives for each value of where Term 1 simplifies to conditional average rates through a shift of indices in the summation where the penultimate step follows because the missing terms are all equal to zero for systems with non-negative abundances: First, consider those reactions for which the step-size of X j is positive, i.e., d kj > 0 then the missing terms correspond to states X j = x j where x j = −d kj , −d jk + 1, . . . , −1. But for all those states the probability to find the system in X j = x j must be equal to zero because they all correspond to negative abundances of X j . Similarly, for reactions for which d jk < 0 the extra terms correspond to states X j = x j where x j = 0, 1, . . . − d jk − 1,. But in all those states the rate of this reaction must be equal to zero because the corresponding transition would take the system to negative abundances of X j .
Considering the above relations at stationarity we thus obtain where the conditional averages are now taken over the stationary state distribution P ss (x), which corresponds to Eq. (2) presented in the main text, and previously derived in (3). Overall, the above derivation is a trivial summation -the simplicity of the argument is obscured through the many indices and keeping track of the variables in general. In any specific system, the summation is readily performed similar to how the moment evolution equations can be readily derived from the master equation using a simple summation over states and shifting of the indices.

Derivation of probability balance relations for specific systems
Our starting point is the probability balance equation which must hold for any given molecular component X i of a complex network that reaches a stationary state in which its probability distribution P ss (x i ) no longer changes. This equations follows trivially from the master equation (see above) and was previously discussed elsewhere (3).
Here, we denoted the states using the index variable m to avoid any confusion with the component itself.
General systems with linear degradation of X 3 . First we derive the invariant balance equation for the following class of systems introduced in the main text.
[S.5] Applying Eq. S.4 to the specific dynamics of Eq. S.5 leads to [S.6] where we have dropped the subscript to denote the stationary state probability distribution for notational convenience. For biochemical reaction systems, for the lowest state m = 0 the relationship simplifies to because molecular numbers cannot be negative (which implies that P (x 3 = −1) = 0). Combining Eq. S.6 with Eq. S.7 then gives rise to the following recursion relation which is Eq. (4) presented in the Main Text.
To get the matrix form of these relations which we use in our optimization algorithm, we re-wrote the conditional averages in Eq. S.8 where P (x 2 , x 3 ) is the joint probability of X 2 and X 3 .
In experiments and computational simulations there will always be a maximum value m max for which P (x 3 = m) = 0 ∀m > m max such that Eq. S.9 can be written in matrix form

Timon Wittenstein, Nava Leibovich, and Andreas Hilfinger
General systems with non-linear degradation of X 3 . The derivation of the probability flux balance relations for systems in which X 3 is degraded as a dimer is equivalent to the one in the previous section. First, we apply Eq. S.4 the specified dynamics of Eq. S.11 to obtain This is a second order recurrence equation that simplifies in the states m = 0 and m = 1 for which the degradation rate is zero, such that for m = 0 Combining these, then leads to the following general balance relation as stated in the main text.

Heuristic approach to determine the regularization parameter
Following the assumption that biochemical reaction rates are reasonably smooth we introduced a regularization term in our optimization in order to penalize large jumps in the reaction rate f ( where the regularization term with Γ ij = δ i,j − 2δ i,j+1 + δ i,j+2 essentially corresponds the discrete second derivative of f (x 2 ). Throughout the paper we used = 1/ √ N which appropriately decreases the importance of the regularization term as the data become more precise. However, it is possible that for other systems another value for the regularization parameter may be advantageous. For completeness, we describe a general method to determine an "optimal" value for for arbitrary cases.
The goal is to pick a regularization parameter which is strong enough to smoothen our estimated reaction rate without introducing too strong a linear bias to the optimization. This can be achieved by determining the value of the unbiased objective function of Eq. S.16 after optimization as a function of the regularization parameter. In Fig. S1A we illustrate the result for the "bistable system" introduced earlier. We see that after a certain value for there is a sharp increase of the objective function which means that the optimization gets heavily biased by the regularization term beyond a certain value. As a heuristic we thus propose to pick a value for the regularization parameter just before the sharp increase. Large values of essentially force the resulting f (  S1. The optimal value for the regularization parameter can be heuristically determined through its effect on the objective function. A) Choosing a regularization parameter strong enough to smoothen the estimated reaction rate without introducing a significant bias leads to the best inference. In this case it was determined to be = 0.01 (orange arrow). B) The regularization is essential to achieve optimal results. Shown here are the inferred rate with the chosen = 0.01 (orange crosses) and without regularization, i.e., = 0 (red crosses). We consider the "noise enhancing" system defined above with changing lifetimes ratios. For such systems, inferring f3(x2) is straightforward when τ2 τ3 and the variability of X2 is slow enough for X3 to adjust to X2-levels and the conditional average x3|x2 directly identifies the production rate of X3. In contrast, when upstream variability in X2 is fast compared to the life-time of X3 the conditional average no longer follows the production rate.  Fig. 4 in the main text we change the shape of the production rate f (x2) = λx n 2 /(x n 2 + K n ) by varying K and n. Shown here is the error of the inferred f (x2) when applying our method with an additional monotonicity constraint. The quality of the inference is significantly improved in almost all regimes.