A multispectral bioluminescence optical tomography algorithm makes use of a partial differential equation (PDE) constrained approach. A sequential quadratic programming (SQP) method is demonstrated that allows for solving both forward and inverse problems at once by updating the forward and inverse variables simultaneously at each step of the optimization iterations. Light propagation in biological tissue is modeled by using the equation of radiative transfer (ERT) and performance of the ERT-based PDE-constrained approach is modeled through numerical and experimental studies.