Due to the small number of copies of molecular species involved, such as DNA, mRNA and regulatory proteins, gene expression is a stochastic phenomenon. In eukaryotic cells, the stochastic effects primarily originate in regulation of gene activity. Transcription can be initiated by a single transcription factor binding to a specific regulatory site in the target gene. Stochasticity of transcription factor binding and dissociation is then amplified by transcription and translation, since target gene activation results in a burst of mRNA molecules, and each mRNA copy serves as a template for translating numerous protein molecules. In the present paper, we explore a mathematical approach to stochastic modeling. In this approach, the ordinary differential equations with a stochastic component for mRNA and protein levels in a single cells yield a system of first-order partial differential equations (PDEs) for two-dimensional probability density functions (pdf). We consider the following examples: Regulation of a single auto-repressing gene, and regulation of a system of two mutual repressors and of an activator-repressor system. The resulting PDEs are approximated by a system of many ordinary equations, which are then numerically solved.