Stochasticity plays an essential role in biochemical systems. Stochastic
behaviors of bimodality, excitability and fluctuations have been observed
in biochemical reaction networks at low molecular numbers. Stochastic
dynamics can be captured by modeling the system using a forward Kolmogorov
equation known in the biochemical literature as the chemical master
equation. The chemical master equation describes the time evolution of the
probability distributions of the molecule species. We develop a stochastic
framework for the design of these time-evolving probability distributions
that includes specifying their uni-/multi-modality, their first moments,
and their rate of convergence to the stationary distribution. By solving
the corresponding optimization programs, we determine the reaction rates
of the biochemical system that satisfy our design specifications. We then
apply the design framework to examples of biochemical reaction networks to
illustrate its strengths and weaknesses.