Many biological networks exhibit stochasticity due to a combinatorial effect of low molecular concentrations and slow system dynamics. One important biological context where stochastic events likely have a large impact is the Bone Morphogenetic Protein (BMP) signaling pathway. BMPs make up the largest subfamily of the Transforming Growth Factor-β
superfamily and are involved in numerous processes including growth, differentiation and diseases [1
]. Due to their potency at driving development, they are also of great value for stem-cell differentiations in cell culture. BMPs activate near maximal signaling at 1nM
concentration, have very slow binding kinetics and require oligomerization between multiple receptor subunits [1
]. These properties naturally lead to conditions for significant and long-duration stochastic fluctuations in cellular signaling. Interestingly, variability of BMP signaling appears to be very low in vivo
, while it is very high in stem cell culture studies [2
]. To understand the differences between in vivo
and in vitro
signaling and determine how various receptor oligomerization events might alter the signal and noise, a more efficient means of solving the steady state distributions for stochastic model was needed that would allow for continuation of both parameters and levels of the BMP pathway components.
Stochastic regulation can negatively impact the robustness of the system [3
] or instead, constructively contribute to the phenotypic variation [5
] in a species. In stochastic reaction networks, the state of a species traverses different trajectories in a probabilistic manner and the distributions of states can be difficult to predict. As more biological data is available, stochastic modeling is becoming increasingly popular to estimate properties in networks where the time evolution of the system is unpredictable and dependent on unavoidable randomness inherent to the system. The complete solution can be calculated from a Chemical Master Equation (CME) [8
], that is based on a Markovian approach that captures the inherent randomness of biochemical systems.
The Chemical Master Equation (CME) describes the dynamics of the probability distribution of a species of chemical reactions. Precisely, the CME captures the rate of change of probability that a system will be in state X
at time t
for all the species of the system. Solution of the CME is practically intractable due to the curse of dimensionality, as the state-space of the system becomes enormously large with increases in the species number and concentrations (number of states nN
, for N
→ species, n
→ copies of each species). Moreover, the system often involves interactions between different time-scales (slow and fast reactions, frequent and infrequent transitions between states) [11
], which add further complexity. Instead, numerical approaches are commonly used [12
] to realize the CME of the stochastic system.
In the analysis of stochastic biochemical networks, steady state probability distributions for each species in the system are determined to measure variability about the deterministic steady state. The deviation around the solution contributes to stochastic noise that can be quantified by measuring the coefficient of variation
(the ratio between the standard deviation σ
and the mean level of species concentration μ
) obtained by solving the CME [9
Frequently, Monte-carlo based simulation approaches [9
] are used to solve stochastic problems. But, there are drawbacks in this approach for networks where the dynamics of different intermediate states of the system are unknown and continuation of several parameters is necessary to explore the system's dependency. As a screen of parameter values becomes necessary in such a scenario, the Monte-carlo based approach doesn't prove to be efficient, as it generally takes longer time to numerically simulate the process and satisfy the imposed conditions. Moreover, simulation times increase with increases in the total number of molecules, species and the number of interactions between species.
In order to ameliorate the computational cost and complexity, we present a method here to approximate the steady state probability distribution by 1) reducing the system's state-space to a finite dimension using truncated state-space method [15
] and 2) subsequently, translating an eigenvalue problem associated with a CME to a system of linear equations. We illustrate that the eigenvector (for an eigenvalue = 0) that represents the steady state probability distribution can be solved algebraically by approximating it as a system of linear equations. Previously, the influence of stochastic fluctuations on system behavior was studied also in [16
], where a moment closure approach was applied to reduce the complexity associated with the identification of distributions. In contrast to the previous studies, here we use a truncated state-space for steady state probability distribution approximation, which is arguably more general since we make no assumptions on the relationships of the moments of the distribution.
The usefulness of the proposed method is illustrated by considering the example networks of BMP signaling, as described earlier in [17
]. Here we examine two potential mechanisms of BMP signaling: 1) regulation between the type I and type II receptors, and 2) regulation by secreted co-factors, such as Crossveinless-2 (Cv-2). First, we apply the approach to the recruitment of Type II receptor into a BMP-bound type I receptor complex to see if such step of receptor oligomerization contributes qualitatively to the noise profile of the system. Following this, we extend our earlier work that focused on extracellular regulation of BMP signaling by SBPs to evaluate the calculation approach and compare results to the Type I/Type II receptor system.
Unlike SBPs, which tend to improve the signal to noise ratio, we did not see a significant change in stochastic variability with inclusion of Type II receptors. This result supports a previous assumption made in [4
] where the recruitment of Type II receptor was excluded to simplify the modeling steps while characterizing the noise profile of a SBP regulated BMP signaling system.