Increasing reliance on complex physiologically-based pharmacokinetic (PBPK) models instead of clinical trials for decision making in drug development calls for an improved framework able to represent a priori parameter dependencies when creating representative virtual populations and performing related statistical inference. We describe here a graph-based solution to harness such dependencies and apply it to inference on hierarchical (population) model parameters. It can model complex parameter and data dependencies, perform efficient Monte Carlo simulations to generate virtual individuals, but also rigorous Bayesian inference on parameters when observed individual covariates are available. Plasma PK data of theophylline were used as an example of application of a PBPK model with built-in covariate structure. A range of models with increasing sophistication and accuracy at describing the data generation process was considered. A stationary MCMC sampler is also described, which has lower complexity than full Bayesian inference. The use of such a sampler and the changes in inference at each stage of the process are discussed.