This paper focuses on the numerical approximation of random lattice reversible Selkov systems. It establishes the existence of numerical invariant measures for random models with nonlinear noise, using the backward Euler-Maruyama (BEM) scheme for time discretization. The study examines both infinite dimensional discrete random models and their corresponding finite dimensional truncations. A classical path convergence technique is employed to demonstrate the convergence of the invariant measures of the BEM scheme to those of the random lattice reversible Selkov systems. As the discrete time step size approaches zero, the invariant measure of the random lattice reversible Selkov systems can be approximated by the numerical invariant measure of the finite dimensional truncated systems.