Bayesian Network (BN) modeling is a prominent and increasingly popular computational systems biology method.It aims at constructing probabilistic networks from the big heterogeneous biological data that reflect the underlying networks of biological relationships. Currently, there is a variety of strategies for evaluating the BN methodology performance, ranging from utilizing the artificial benchmark datasets and models, to the specialized biological benchmark datasets, to the simulation studies that generate synthetic data from the predefined network models. The latter is arguably the most comprehensive approach; however, existing implementations are typically limited by their reliance on the SEM (structural equation modeling) framework, with many explicit and implicit assumptions that might not be realistic in a typical biological data analysis scenario. In this study, we develop an alternative, purely probabilistic, simulation framework that is a more appropriate fit with the real biological data and biological network models. In conjunction, we also expand on our current understanding of the theoretical notions of causality and dependence / conditional independence in BNs and Markov Blankets within.