Abstract Mathematical models of biomolecular networks are commonly used to study cellular processes; however, their usefulness to explain and predict dynamic behaviors is often questioned due to the unclear relationship between parameter uncertainty and network dynamics. In this work, we introduce PyDyNo (Python Dynamic analysis of biochemical NetwOrks), a non-equilibrium reaction-flux based analysis to identify dominant reaction paths within a biochemical reaction network calibrated to experimental data. We first show, in a simplified apoptosis execution model, that Bayesian parameter optimization can yield thousands of parameter vectors with equally good fits to experimental data. Our analysis however enables us to identify the dynamic differences between these parameter sets and identify three dominant execution modes. We further demonstrate that parameter vectors from each execution mode exhibit varying sensitivity to perturbations. We then apply our methodology to JAK2/STAT5 network in colony-forming unit-erythroid (CFU-E) cells to identify its signal execution modes. Our analysis identifies a previously unrecognized mechanistic explanation for the survival responses of the CFU-E cell population that would have been impossible to deduce with traditional protein-concentration based analyses. Impact Statement Given the mechanistic models of network-driven cellular processes and the associated parameter uncertainty, we present a framework that can identify dominant reaction paths that could in turn lead to unique signal execution modes (i.e., dominant paths of flux propagation), providing a novel statistical and mechanistic insights to explain and predict signal processing and execution.

Emergent signal execution modes in biochemical reaction networks calibrated to experimental data