Spatiotemporal control of the complex chemical microenvironment is of great importance to many fields within life science. One way to facilitate such control is to construct delivery circuits, comprising arrays of dispensing outlets, for ions and charged biomolecules based on ionic transistors. This allows for addressability of ionic signals, which opens up for spatiotemporally controlled delivery in a highly complex manner. One class of ionic transistors, the ion bipolar junction transistors (IBJTs), is especially attractive for these applications because these transistors are functional at physiological conditions and have been employed to modulate the delivery of neurotransmitters to regulate signaling in neuronal cells. Further, the first integrated complementary ionic circuits were recently developed on the basis of these ionic transistors. However, a detailed understanding of the device physics of these transistors is still lacking and hampers further development of components and circuits. Here, we report on the modeling of IBJTs using Poisson's and Nernst-Planck equations and the finite element method. A two-dimensional model of the device is employed that successfully reproduces the main characteristics of the measurement data. On the basis of the detailed concentration and potential profiles provided by the model, the different modes of operation of the transistor are analyzed as well as the transitions between the different modes. The model correctly predicts the measured threshold voltage, which is explained in terms of membrane potentials. All in all, the results provide the basis for a detailed understanding of IBJT operation. This new knowledge is employed to discuss potential improvements of ion bipolar junction transistors in terms of miniaturization and device parameters.