Pharmacokinetic/Pharmacodynamic (PKPD) models are empirical models of the physiological and pharmacological systems often used to describe the kinetics and behavior of drugs in the human body. Nonlinear Mixed Effects (NLME) statistical methods help identify the parameters of the PKPD models and quantify the differences between individuals by integrating models at the population and individual scales. In this talk I introduce PuMaS.jl, a Julia based software for simulating and estimating PKPD, physiology based PK (PBPK), quantitative systems pharmacology (QSP), etc. models used in pharmacology. I will begin by describing approximations to the marginal likelihood which are used to make the quantities efficiently computable and demonstrate on real data how these models can be fit with Optim.jl to reveal population-level characteristics. Additionally, I will demonstrate the ability to utilize DynamicHMC.jl to perform Bayesian estimation of population and individual parameters. Together, this demonstrates a Julia-based data-driven approach to handle complex problems in individualizing dosing.