This paper describes a numerical approach to the modeling of PNP HBTs in the InP-based materials systems (InP/InGaAs and InAlAs/InGaAs). Initial device analysis was achieved in the drift-diffusion limit by self-consistent numerical solution of the Poisson, carrier continuity and conductor equations subject to the device's geometry and boundary conditions imposed by the device's biasing. Simulation results are compared with the available experimental results and good agreement is found. For the (\operatorname{InP} / \operatorname{InGaAs}) and (\mathrm{InAlAs} / \mathrm{InGaAs}) heterojunctions, the valence band discontinuities are larger than for the (\mathrm{AlGaAs} / \mathrm{GaAs}) system so grading of the emitter-base junction and tunneling effects are important. For completeness, nonclassical effects were also considered. For the emitter-base junctions, hole tunneling was considered, particularly at low forward bias. The inverse dependence of the hole tunneling on effective mass was found to lead to significantly more light hole than heavy hole tunneling in calculating the emitter injection current. In addition, since very narrow base regions ( (25-35 \mathrm{~nm}) ) can be employed while keeping the base spreading resistance low due to the electron's higher mobility, ballistic hole transport should also be considered.