Background In one way or the other, all modern parametrizations of the nuclear energy density functional (EDF) do not respect the exchange symmetry associated with Pauli's principle. It has been recently shown that this practice jeopardizes multi-reference (MR) EDF calculations by contaminating the energy with spurious self-interactions that, for example, lead to finite steps or even divergences when plotting it as a function of collective coordinates [J. Dobaczewski et al., Phys. Rev. C 76, 054315 (2007); D. Lacroix et al., Phys. Rev. C 79, 044318 (2009)]. As of today, the only viable option to bypass these pathologies is to rely on EDF kernels that enforce Pauli's principle from the outset by strictly and exactly deriving from a genuine, i.e. density-independent, Hamilton operator.Purpose The objective is to build cutting-edge parametrizations of the EDF kernel deriving from a pseudo potential that can be safely employed in symmetry restoration and configuration mixing calculations.Methods We wish to develop the most general Skyrme-like EDF parametrization containing linear, bilinear and trilinear terms in the density matrices with up to two gradients, under the key constraint that it derives strictly from an effective Hamilton operator. While linear and bilinear terms are obtained from a standard one-body kinetic energy operator and a (density-independent) two-body Skyrme pseudo-potential, the most general three-body Skyrme-like pseudo-potential containing up to two gradient operators is constructed to generate the trilinear part. The present study is limited to central terms. Spin-orbit and tensor will be addressed in a forthcoming paper.Results The most general central Skyrme-type zero-range three-body interaction is built up to second order in derivatives. The complete trilinear energy density functional, including time-odd and T = 1 pairing parts, is derived along with the corresponding normal and anomalous fields entering the Hartree-Fock-Bogoliubov equations of motion. Its building blocks are the same local densities that the standard Skyrme functional is constructed from. The central three-body pseudopotential is defined out of six independent parameters. Expressions for bulk properties of symmetric, isospin-asymmetric and spin-polarized homogeneous nuclear matter, as well as associated Landau parameters, are given.Conclusions This study establishes a first step towards a new generation of nuclear energy density functionals that respect Pauli's principle and that can be safely used in predictive and spuriousity-free SR and MR EDF calculations.