We use the `Evolution and assembly of galaxies and their environments' (EAGLE) cosmological simulation to investigate the effect of baryons on the density profiles of rich galaxy clusters. We focus on EAGLE clusters with M200 > 1014 M⊙ of which we have six examples. The central brightest cluster galaxies (BCGs) in the simulation have steep stellar density profiles, ρ*(r) ∝ r-3. Stars dominate the mass density for r < 10 kpc, and, as a result, the total mass density profiles are steeper than the Navarro-Frenk-White (NFW) profile, in remarkable agreement with observations. The dark matter halo itself closely follows the NFW form at all resolved radii (r ≳ 3.0 kpc). The EAGLE BCGs have similar surface brightness and line-of-sight velocity dispersion profiles as the BCGs in the sample of Newman et al., which have the most detailed measurements currently available. After subtracting the contribution of the stars to the central density, Newman et al. infer significantly shallower slopes than the NFW value, in contradiction with the EAGLE results. We discuss possible reasons for this discrepancy, and conclude that an inconsistency between the kinematical model adopted by Newman et al. for their BCGs, which assumes isotropic stellar orbits, and the kinematical structure of the EAGLE BCGs, in which the orbital stellar anisotropy varies with radius and tends to be radially biased, could explain at least part of the discrepancy.