The emerging field of Glycomics requires the development of systems-based modeling strategies to relate glycosyltransferase gene expression and enzyme activity with carbohydrate structure and function. We describe the application of object oriented programming concepts to define glycans, enzymes, reactions, pathways and compartments for modeling cellular glycosylation reaction networks. These class definitions are combined with current biochemical knowledge to define potential reaction networks that participate in the formation of the sialyl Lewis-X (sLe(X)) epitope on O-glycans linked to a leukocyte cell-surface glycoprotein, P-selectin Glycoprotein Ligand-1 (PSGL-1). Subset modeling, hierarchical clustering, principal component analysis and adjoint sensitivity analysis are applied to refine the reaction network and to quantify individual glycosyltransferase rate constants. Wet-lab experiments validate estimates from computer modeling. Such analysis predicts that sLe(X) expression varies directly with sialyltransferase alpha2,3ST3Gal-IV expression and inversely with alpha2,3ST3Gal-I/II. SBML files for all converged models are available at http://www.eng.buffalo.edu/~neel/bio_reaction_network.html