We construct nuclear energy density functionals in terms of derivatives of densities up to sixth order. A phenomenological functional built in this way conforms to the ideas of the density matrix expansion and is rooted in the expansions characteristic to effective theories. It builds on the standard functionals related to the contact and Skyrme forces, which constitute the zero-order and second-order expansions, respectively. At 6th order, the full functional with density-independent coupling constants, and with the isospin degree of freedom taken into account, contains 376 terms. Imposing different symmetries gives relations between the terms and reduces the number of free parameters to reasonable values.