We use molecular dynamics simulations with a dissipative particle dynamics (DPD) thermostat to study the behavior of nanosized inclusions (colloids) in a polymer brush which is in contact with an explicit solvent in the NPT ensemble. The brush is described by a bead-spring model for flexible polymer chains, grafted on a solid substrate, while the polymer-soluble nanoparticles in the solution are taken as hard spheres. By varying the chain length N, the grafting density of the brush, σ g , and the size of the nanoparticles b, we determine the equilibrium particle penetration depth δ and the average concentration of nanoinclusions ϕ ¯ nano in the penetration layer δ at constant pressure. In agreement with a recent theoretical prediction, we demonstrate that for nanoinclusions of size b ⩾ b ∗ the thickness of this layer δ ∝ h ( b ∗ / b ) 3 where h is brush height and b ∗ ∝ σ g - 2 / 3 is a typical size below which smaller particles are uniformly distributed in the brush. We also observe that particles, larger than some threshold value b max do not mix with the brush. The mean density of nanoinclusions is found to scale as ϕ ¯ nano ∝ ( b ∗ / b ) 3 within the whole range of parameter variation. The diffusivity of nanoparticles, embedded in the polymer brush, in direction perpendicular to the grafting plane is found to be up to 20% higher than parallel to the plane. The variation of the respective diffusion coefficients D ⊥ nano and D | | nano changes with growing volume fraction of the nanoparticles in agreement with theoretical predictions.