In clinical trials, the response of a given subject often depends on the selected treatment as well as on some covariates. We study optimal approximate designs of experiments in the models with treatment and covariate effects. We allow for the variances of the responses to depend on the chosen treatments, which introduces heteroscedasticity into the models. For estimating systems of treatment contrasts and linear functions of the covariates, we extend known results on D-optimality of product designs by providing product designs that are optimal with respect to general eigenvalue-based criteria. In particular, A- and E-optimal product designs are obtained. We then formulate a method based on linear programming for constructing optimal designs with smaller supports from the optimal product designs. The sparser designs can be more easily converted to practically applicable exact designs. The provided results and the proposed sparsification method are demonstrated on some examples.