Function on scalar regression models relate functional outcomes to scalar predictors through the conditional mean function. With few and limited exceptions, many functional regression frameworks operate under the assumption that covariate information does not affect patterns of covariation. In this manuscript, we address this disparity by developing a Bayesian functional regression model, providing joint inference for both the conditional mean and covariance functions. Our work hinges on basis expansions of both the functional evaluation domain and covariate space, to define flexible non-parametric forms of dependence. To aid interpretation, we develop novel low-dimensional summaries, which indicate the degree of covariate-dependent heteroskedasticity. The proposed modeling framework is motivated and applied to a case study in functional brain imaging through electroencephalography, aiming to elucidate potential differentiation in the neural development of children with autism spectrum disorder.