Investigation of HIV viral dynamics is important for understanding the HIV pathogenesis and for development of treatment strategies. Perelson et al. demonstrated that simple viral dynamic models fit to data on viral load as measured by plasma HIV-RNA could produce estimates of rates of clearance of virus and of infected CD4+ T-lymphocytes. In this paper we extend the work of Perelson et al. by proposing models with less restrictive assumptions about drug activity. Our models take into account the fact that infectious and non-infectious virions are produced by infected T-cells both before and after the treatment. We also show that direct measurement of infectious virus load provides sufficient information for estimation of antiretroviral drug efficacy parameter. For characterizing viral dynamics of populations and estimation of dynamic parameters, we propose a hierarchical non-linear model. Compared to other methods such as the non-linear least square method used by Perelson et al., we show that the proposed approach has the following advantages: (i) it is more appropriate for modelling within-patient and between-patient variation and to characterize the population dynamics; (ii) it is flexible enough to deal with both rich and sparse individual data; (iii) it has more power to detect model misspecification; (iv) it allows incorporation of covariates for viral dynamic parameters; (v) it makes more efficient use of between-subject information to get better parameter estimates. We give two simulation examples to illustrate the proposed approach and its advantages. Finally, we discuss practical issues regarding the clinical trial design for viral dynamic studies.