Gels are comprised of polymer networks swelled by some interstitial solvent. They are under wide investigation by material scientists and engineers for their broad applicability in fields ranging from adhesives to tissue engineering. Gels’ mechanical properties greatly influence their efficacy in such applications and are largely dictated by their underlying microstructures and constituent-scale properties. Yet predictively mapping the local-to-global property functions of gels remains difficult due - in part - to the complexity introduced by solute-solvent interactions. We here introduce a novel, discrete mesoscale modeling method that preserves local solute concentration-dependent gradients in osmotic pressure through the Flory-Huggins mixing parameter, χ. The iteration of the model used here replicates gels fabricated from telechelically crosslinked star-shaped polymers and intakes χ, macromer molecular weight (Mw), crosslink functionality (f), and as-prepared solute concentration (ϕ*) as its inputs, all of which are analogues to the control parameters of experimentalists. Here we demonstrate how this method captures solvent-dependent homogenization (χ≤0.5) or phase separation (χ>0.5) of polymer suspensions in the absence of phenomenological pairwise potentials. We then demonstrate its accurate, ab initio prediction of gel topology, isotropic swelling mechanics, and uniaxial tensile stress for a 10k tetra-PEG gel. Finally, we use the model to predict trends in the mechanical response and failure of multi-functional PEG-based gels over a range of Mw and f, while investigating said trends’ micromechanical origins. The model predicts that increased crosslink functionality results in higher initial chain stretch (as measured at the equilibrated swollen state) for gels of the same underlying chain length, which improves modulus and failure stress but decreases failure strain and toughness.