The idea of describing a many-electron system using only its electron density, i.e. without constructing its wavefunction, was initiated in the works of Thomas and Fermi. Hohenberg-Kohn theorems of modern density functional theory transformed this idea into an exact theory. The Kohn-Sham formalism, widely used in computer simulations of polyatomic systems today, is based on these theorems but is not orbital-free. It reintroduces orbitals to minimize errors in approximating the total energy. The present review concerns an alternative formalism based also on Hohenberg-Kohn theorems, in which orthogonal orbitals are used not for the whole system but only for subsystems [Cortona, Phys. Rev. B, 44 (1991) 8454]. These orbitals are derived from Kohn-Sham-like one-electron equations, called here Kohn-Sham Equations with Constrained Electron Density (KSCED), in which all terms representing the interactions between the subsystems are expressed as universal functionals of electron density. This formulation provides the formal basis for the orbital-free embedding in first-principles based multi-level simulations of complex systems, in which the orbital-level is retained for a selected subsystem, whereas its environment is described at the orbital-free level [Wesolowski and Warshel, J. Phys. Chem., 97 (1993) 8050]. The formal aspects, development and testing of relevant approximate density functionals, and the possible use of the orbital-free embedding in multi-level modelling are covered in detail in this review. Examples of applications, especially those concerning the electronic structure of embedded systems in the condensed phase are provided.