Polymerization and formation of cross-linked polymer networks are important processes in manufacturing, materials fabrication, and in the case of hydrated polymer networks, synthesis of biomedical materials, drug delivery, and tissue engineering. 'While considerable research has been devoted to the modeling of polymer networks to determine averaged, mean-field, global properties, there are fewer studies that specifically examine the variance of the composition across "microregions" (composed of a large, but finite, number of polymer network strands) within the larger polymer network. Here, we mathematically model the stochastic formation of polymer networks comprised of linear network strands with two identical reactive end-groups that undergo an endlinking gelation process. We focus on networks formed by chain-growth polymerization but also give examples of how our model can be extended to networks formed by step-growth polymerization. We introduce a master equation that describes the evolution of the probabilities of possible network microregion configurations as a function of time and the extent of reaction. We specifically focus on the dynamics of network formation and the statistical variability of the gel microregions, particularly at intermediate extents of reaction. We also consider possible annealing effects and study how cooperative binding between the two end-groups on a single network strand affects network formation. Our results allow for a more detailed and thorough understanding of polymer network dynamics and variability of network properties.