Individuals differ in how they respond to a given treatment, and characterizing such variability in treatment response is an imperative aim in personalized medicine. We propose a general method to model treatment response heterogeneity, through identification of treatment-covariate interactions honoring different hierarchy conditions. We construct a single-step $l_1$ norm penalty procedure that maintains the hierarchical structure of interactions in a sense that the treatment-covariate interaction term is included in the model only when either of the covariate or the covariate and the treatment both have non-zero main effects. We explore several parameterization schemes with different constraints added to Lasso that enforce the hierarchical interaction restriction. We solve the resulting constrained optimization problem using a spectral projected gradient method. We compare our methods to the unstructured Lasso using simulation studies covering a variety of scenarios for treatment-covariate interactions. The simulations show that our methods yield more parsimonious models and outperform unstructured Lasso in terms of prediction performance, and in terms of the ability to correctly identify non-zero treatment covariate interactions. The superior performance of our methods are also corroborated by an application on a large randomized clinical trial data investigating a drug for treating congestive heart failure (N=2,569). Our methods can be applied to continuous, binary and time to event outcome, providing a well-suited approach with sufficient flexibility in terms of parameterization for doing secondary analysis in clinical trials to analyze heterogeneity treatment effect.