  Reaction networks are systems in which the populations of a finite number of
species evolve through predefined interactions. Such networks are found as
modeling tools in many biological disciplines such as biochemistry, ecology,
epidemiology, immunology, systems biology and synthetic biology. It is now
well-established that, for small population sizes, stochastic models for
biochemical reaction networks are necessary to capture randomness in the
interactions. The tools for analyzing such models, however, still lag far
behind their deterministic counterparts. In this paper, we bridge this gap by
developing a constructive framework for examining the long-term behavior and
stability properties of the reaction dynamics in a stochastic setting. In
particular, we address the problems of determining ergodicity of the reaction
dynamics, which is analogous to having a globally attracting fixed point for
deterministic dynamics. We also examine when the statistical moments of the
underlying process remain bounded with time and when they converge to their
steady state values. The framework we develop relies on a blend of ideas from
probability theory, linear algebra and optimization theory. We demonstrate that
the stability properties of a wide class of biological networks can be assessed
from our sufficient theoretical conditions that can be recast as efficient and
scalable linear programs, well-known for their tractability. It is notably
shown that the computational complexity is often linear in the number of
species. We illustrate the validity, the efficiency and the wide applicability
of our results on several reaction networks arising in biochemistry, systems
biology, epidemiology and ecology. The biological implications of the results
as well as an example of a non-ergodic biological network are also discussed.
