Analytical expressions for the saturation density of asymmetric nuclear matter as well as its binding energy and incompressibility at saturation density are given up to fourth order in the isospin asymmetry δ = (ρ n − ρ p )/ρ using 11 characteristic parameters defined by the density derivatives of the binding energy per nucleon of symmetric nuclear matter, the symmetry energy E sym (ρ), and the fourth-order symmetry energy E sym,4 (ρ) at normal nuclear density ρ 0 . Using an isospin-and momentum-dependent modified Gogny interaction (MDI) and the Skyrme-Hartree-Fock (SHF) approach with 63 popular Skyrme interactions, we have systematically studied the isospin dependence of the saturation properties of asymmetric nuclear matter, particularly the incompressibilityat saturation density. Our results show that the magnitude of the higher order K sat,4 parameter is generally small compared to that of the K sat,2 parameter. The latter essentially characterizes the isospin dependence of the incompressibility at saturation density and can be expressed asL, where L and K sym represent, respectively, the slope and curvature parameters of the symmetry energy at ρ 0 and J 0 is the third-order derivative parameter of symmetric nuclear matter at ρ 0 . Furthermore, we have constructed a phenomenological modified Skyrme-like (MSL) model that can reasonably describe the general properties of symmetric nuclear matter and the symmetry energy predicted by both the MDI model and the SHF approach. The results indicate that the higher order J 0 contribution to K sat,2 generally cannot be neglected. In addition, it is found that there exists a nicely linear correlation between K sym and L as well as between J 0 /K 0 and K 0 . These correlations together with the empirical constraints on K 0 , L, E sym (ρ 0 ), and the nucleon effective mass lead to an estimate of K sat,2 = −370 ± 120 MeV.