Glutathione reductase (GR) catalyzes the reduction of oxidized glutathione (GSSG) to reduced glutathione (GSH) using NADPH as the reducing cofactor, and thereby maintains a constant GSH level in the system. GSH scavenges superoxide (O2·−) and hydroxyl radicals (OH·) non-enzymatically or by serving as an electron donor to several enzymes involved in reactive oxygen species (ROS) detoxification. In either case, GSH oxidizes to GSSG and is subsequently regenerated with the catalytic action of GR. Though GR kinetic mechanism has been extensively studied under different experimental conditions with variable substrates and products, the catalytic mechanism has not been studied in terms of a mechanistic model that accounts for the effects of the substrates and products on the reaction kinetics. The aim of the current study is therefore to develop a comprehensive mathematical model for the catalytic mechanism of GR. We use available experimental data on GR kinetics from various species/sources to develop the mathematical model and estimate the associated model parameters. The model simulations are consistent with the experimental observations that GR operates via both ping-pong and sequential branching mechanisms based on relevant concentrations of its reaction substrate GSSG. Furthermore, we show the observed pH-dependent substrate-inhibition of GR activity by GSSG and bi-modal behavior of GR activity with pH. The model presents a unique opportunity to understand the effects of products on the kinetics of GR. The model simulations show that under physiological conditions, where both substrates and products are present, the flux distribution depends on the concentrations of both GSSG and NADP+ with ping-pong flux operating at low levels and sequential flux dominating at higher levels. The kinetic model of GR may serve as a key module for the development of integrated models for ROS scavenging system to understand protection of cells under normal and oxidative stress conditions.