Recent experiments have unambiguously established that biological systems can have significant cell to cell variations in gene expression levels even in isogenic populations. Computational approaches to studying gene expression in cellular systems should capture such biological variations for a more realistic representation. We present a simple, fully probabilistic approach to the modeling of gene regulatory networks that allows for fluctuations in the gene expression levels. The new algorithm uses a very simple representation for the genes, and accounts for the repression or induction of the genes and for the biological variations among isogenic populations simultaneously. We have tested the new algorithm on a recently bioengineered synthetic gene network library, and found a good agreement between model predictions and experimental results.