We aim to develop a new model reduction approach for biochemical reaction networks (BCRN) governed by mass action kinetics using data-driven techniques. In this framework, we first develop a compact discrete mathematical formulation for modelling BCRN which is based on the Laplacian of the graph of complexes. This graph has complexes, i.e., the left- and right-hand sides of the reactions as nodes and reactions as edges. The next step in the proposed model reduction procedure is to select a set of species concentrations to be eliminated from the model and to link the remaining species concentrations via possible networks that obey mass action kinetics. We then optimize the difference between the original and the reduced model over all the possible networks using large time series concentration datasets generated from the original model to determine a reduced model which is closest to the original. This optimization problem is solved using neural networks.