Given the optimal interpolation points~$\sigma_1,\dots ,\sigma_r$, it is well-known how to obtain the ℋ₂-optimal reduced order model of order~$r$ for a linear time-invariant system of order $n\gg r$. Our approach to linear time-invariant systems depending on parameters~$p$ is to approximate their parametric dependence as a so-called metamodel, which in turn allows us to set up the corresponding parametrized reduced order models. The construction of the metamodel we suggest involves the coefficients of the characteristic polynomial together with $k$-means clustering and radial basis function interpolation, and thus allows for an accurate and efficient approximation of~$\sigma_1(p),\dots,\sigma_r(p)$. As the computation still includes large system solves, this metamodel is not sufficient to construct a fast and truely parametric reduced system. Setting up a medium size model without extra cost, we present a possible answer to this. We illustrate the proposed method with several numerical examples.