  Stochasticity plays important roles in reaction systems. Vector fields of
probability flux and velocity characterize time-varying and steady-state
properties of these systems, including high probability paths, barriers,
checkpoints among different stable regions, as well as mechanisms of dynamic
switching among them. However, conventional fluxes on continuous space are
ill-defined and are problematic when at boundaries of the state space or when
copy numbers are small. By re-defining the derivative and divergence operators
based on the discrete nature of reactions, we introduce new formulations of
discrete fluxes. Our flux model fully accounts for the discreetness of both the
state space and the jump processes of reactions. The reactional discrete flux
satisfies the continuity equation and describes the behavior of the system
evolving along directions of reactions. The species discrete flux directly
describes the dynamic behavior in the state space of the reactants such as the
transfer of probability mass. With the relationship between these two fluxes
specified, we show how to construct time-evolving and steady-state global
flow-maps of probability flux and velocity in the directions of every species
at every microstate, and how they are related to the outflow and inflow of
probability fluxes when tracing out reaction trajectories. We also describe how
to impose proper conditions enabling exact quantification of flux and velocity
in the boundary regions, without the difficulty of enforcing artificial
reflecting conditions. We illustrate the computation of probability flux and
velocity using three model systems, namely, the birth-death process, the
bistable Schl\"ogl model, and the oscillating Schnakenberg model.
