Hydrogels have seen widespread application across biomedical sciences and there is considerable interest in using hydrogels, including agarose, for creating in vitro three-dimensional environments to grow cells and study mechanobiology and mechanotransduction. Recent advances in the preparation of agarose gels enable successful encapsulation of viable cells at gel concentrations as high as 5%. Agarose with a range of gel concentrations can thus serve as an experimental model mimicking changes in the 3-D microenvironment of cells during disease progression and can facilitate experiments aimed at probing the corresponding mechanobiology, e.g. the evolving mechanobiology of chondrocytes during the progression of osteoarthritis. Importantly, whether stresses (forces) or strains (displacements) drive mechanobiology and mechanotransduction is currently unknown. We can use experiments to quantify mechanical properties of hydrogels, and imaging to estimate microstructure and even strains; however, only computational models can estimate intra-gel stresses in cell-seeded agarose constructs because the required in vitro experiments are currently impossible. Finite element modeling is well-established for (computational) mechanical analyses, but accurate constitutive models for modeling the 3-D mechanical environments of cells within high-stiffness agarose with varying gel concentrations are currently unavailable. In this study we aimed to establish a 3-D constitutive model of high-stiffness agarose with a range of gel concentrations. We applied a multi-step, physics-based optimization approach to separately fit subsets of model parameters and help achieve robust convergence. Our constitutive model, fitted to experimental data on progressive stress-relaxations, was able to predict reaction forces determined from independent experiments on cyclical loading. Our model has broad applications in finite element modeling aimed at interpreting mechanical experiments on agarose specimens seeded with cells, particularly in predicting distributions of intra-gel stresses. Our model and fitted parameters enable more accurate finite element simulations of high-stiffness agarose constructs, and thus better understanding of experiments aimed at mechanobiology, mechanotransduction, or other applications in tissue engineering.