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) nonenzymatically 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 by the catalytic action of GR. Although the GR kinetic mechanism has been extensively studied under various 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 this 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 observation 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 bimodal 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 systems to understand protection of cells under normal and oxidative stress conditions.