Anomalous diffusion is an established phenomenon but still a theoretical challenge in non-equilibrium statistical mechanics. Physical models are built incrementally, and the most recent and most general family is based on the fractional Brownian motion (fBm) with a random diffusion coefficient (superstatistical fBm) together with a time-dependent random Hurst parameter. We provide here a dynamical foundation for such general family of models. We consider a dynamical system describing the motion of a test-particle surrounded by $N$ Brownian particles with different masses. This dynamic is governed by underdamped Langevin equations. Physical principles of conservation of momentum and energy are met. We prove that, in the limit $N\to\infty$, the test-particle diffuses in time according to a quite general (non-Markovian) Gaussian process whose covariance function is determined by the distribution of the masses of the surround-particles. In particular, with proper choices of the distribution of the masses of the surround-particles, we obtain fBm together with a number of other special cases of interest in modelling anomalous diffusion including time-dependent anomalous exponent. Furthermore, when the ensemble heterogeneity of the surround-particles embodying the environment becomes non-uniform and joins with the individual inhomogeneity of the test-particles, we show that, in the limit $N\to\infty$, the test-particle diffuses in time according to a quite general conditionally Gaussian process that can be calibrated into a fBm with random diffusion coefficient and random time-dependent Hurst parameter. We conclude our study by reporting the generalised Kolmogorov--Fokker--Planck equations associated to these highly general processes.