Genistein is an endocrine-active compound found naturally in soy products. It has been linked to various health effects, both beneficial and adverse. The liver is a major site of genistein transformation. Experimental data suggest genistein is metabolized in the liver into its glucuronide form, and taken back into the gut lumen via biliary excretion. The data show a nontrivial, dose dependent delay in biliary excretion of genistein. Traditional physiologically-based pharmacokinetic (PBPK) modeling methods fail to accurately describe the observed data. We have developed several models that incorporate techniques not typically found in PBPK models to simulate the observed dynamics. The first of these models developed is based on delay differential equations (DDE), where the observed lag in biliary excretion is mathematically described by a time delay. Existence and uniqueness of a solution to the system of equations was obtained and the unknown parameters were obtained via an inverse problem formulation. The nonlinear system of state-dependent delay equations was approximated by a three-stage, implicit, Runge-Kutta method. Using a statistical hypothesis test, we showed that the delay model with the optimal set of parameters obtained is a statistically significant improvement over the PBPK model in simulating the experimental data. Our second modeling approach was taken by considering a dispersion modeling technique. We seeked to develop a fully functioning model for the liver that simulates the distribution, metabolism and excretion of chemicals, and is able to accommodate spatial variations in biologically- and physiologically-based parameters. Our modeling strategy considers the liver as a series of cylindrical tubes, one representing the blood vessel space and one representing the bile duct space, with hepatocytes in between. Dispersion coefficients were adjusted to create a biologically relevant distribution of the concentration of genistein and its metabolites. We established existence and uniqueness of a solution to the system of equations and continuous dependence of the solution on the initial data. A numerical code based on the finite element method was developed to solve the mixed system of nonlinear ordinary and partial differential equations. Parameter estimations were obtained via an inverse problem formulation.