|Abstract:|| This work addresses the problem of performing functional linear regression when the response variable is represented as a probability density function (PDF). PDFs are interpreted as functional compositions, that are objects carrying primarily relative information. In this context, the unit integral constraint allows to single out one of the possible representations of a class of equivalent measures. On these bases, a function-on-scalar regression model with distributional response is proposed, by relying on the theory of Bayes Hilbert spaces. The geometry of Bayes spaces allows capturing all the key inherent feature of distributional data (e.g., scale invariance, relative scale). A B-spline basis expansion combined with a functional version of the centred log-ratio transformation is employed for actual computations. For this purpose, a new key result is proved to characterize B-spline representations in Bayes spaces. We show the potential of the methodological developments on a real case study, dealing with metabolomics data. Here, a bootstrap-based study is also performed for the uncertainty quantification of the obtained estimates.