  In this paper, we provide a detailed theoretical analysis of the numerical
scheme introduced in J. Comput. Phys. 436 (2021) 110253 for the reaction
kinetics of a class of chemical reaction networks that satisfies detailed
balance condition. In contrast to conventional numerical approximations, which
are typically constructed based on ordinary differential equations (ODEs) for
the concentrations of all involved species, the scheme is developed using the
equations of reaction trajectories, which can be viewed as a generalized
gradient flow of physically relevant free energy. The unique solvability,
positivity-preserving, and energy-stable properties are proved for the general
case involving multiple reactions, under a mild condition on the stoichiometric
matrix.
