Most biochemical reactions in living cells are open system interacting with environment through chemostats. At a mesoscopic scale, the number of  each species in those biochemical reactions can be modeled by the random time-changed Poisson processes. To characterize the macroscopic behaviors in the large volume limit, the law of large number in path space determines a mean-field limit nonlinear Kurtz ODE, while the WKB expansion yields a Hamilton-Jacobi equation and the corresponding Lagrangian gives the good rate function in the large deviation principle. A parametric variation principle can be formulated to compute the reaction paths. We propose a gauge-symmetry criteria for a class of non-equilibrium chemical reactions including  enzyme reactions, which identifies a new concept of balance within the same reaction vector due to flux grouping degeneracy. With this criteria,  we (i) formulate an Onsager-type gradient flow structure in terms of the energy landscape given by a steady solution to the Hamilton-Jacobi equation;  (ii) find transition paths between multiple non-equilibrium steady states (rare events in biochemical reactions).  We illustrate this idea through a bistable catalysis  reaction. In contrast to the standard diffusion approximations via Kramers-Moyal expansion, a new drift-diffusion approximation sharing the same gauge-symmetry is constructed based on the Onsager-type gradient flow formulation to compute the correct energy barrier.