Truncated series expansions for the species partial molar fugacity coefficients in ternary and quaternary dilute systems are derived for the systematic study of mixed solutes in pure and mixed solvents at near critical conditions.  Then, explicit molecular-based expressions for the expansion coefficients are drawn in terms of direct and total correlation function integrals associated with the actual system's microstructure. Finally, these self-consistent formal expressions are used (a) to interpret the potential synergistic effects on the solubility enhancement caused by the presence of either a co-solute, a co-solvent, or both in terms of the solvation behavior of the dilute species, and (b) to highlight, and illustrate with examples from literature, some frequent pitfalls in the molecular modeling of these mixtures leading to serious thermodynamic inconsistencies.  This research was sponsored by the Division of Chemical Sciences, Geosciences, and Biosciences, Office of Basic Energy Sciences under contract number DE-AC05-00OR22725 with Oak Ridge National Laboratory, managed and operated by UT-Battelle, LLC.