Classic hemodialysis only provides a limited removal of protein bound uremic toxins (PBUT) in patients with chronic kidney disease. A bioartificial kidney device, BAK, composed of a living cell monolayer of conditionally immortalized proximal tubule epithelial kidney cells (ciPTEC) cultured of hollow fiber polymeric membrane can remove protein bound uremic toxins from the blood in combination with classic hemodialysis. The development and clinical implementation of the BAK requires lots of optimization. This investigation is expensive and time consuming therefore modeling studies could help to optimize experiments and improve its design. In this work, a 3D mathematical model of the BAK is developed. The transport and reaction mechanisms associated with the removal of PBUT indoxyl sulfate are considered and various conditions are simulated. The model describes a single hollow fiber membrane and considers different domains for the blood flow, the membrane, the cell monolayer, and the dialysate region. A mathematical description of the relevant transport and/or reaction mechanisms is provided in each domain, and the corresponding differential equations are solved numerically. Since not all the modeling constants are experimentally available, a parametric study is performed for their quantification, including the active transport kinetics of the toxins through the cell monolayer, in comparison to the passive transport rates by diffusion. The parametric study also provides a background for the extraction of usually unknown quantities, including notably the Organic Anion Transporter (OAT) concentrations, with the support of experimental data. Satisfactory reproduction of experimental findings is achieved, and the role of systemic variables that affect significantly the uremic toxin removal is identified.