Optoacoustic (OA) imaging combined with reversibly photoswitchable proteins has emerged as a promising technology for the high-sensitivity and multiplexed imaging of cells in live tissues in preclinical research. Through carefully designed illumination schedules of ON and OFF laser pulses, the resulting OA signal is a multiplex of different reporter species and the background. We propose a model-based variational framework to computationally unmix and image different species of photo-switching reporters using optoacoustic tomography. It is based on a detailed mathematical description of the photo-switching mechanism, which models how relevant physical parameters such as the kinetic constants and light fluence impact the switching signal. We introduce an algorithm that operates on images, as opposed to traditional pixelwise approaches. It takes the form of an iterative inversion combined with tailored ℓ1 and total-variation regularization to increase the robustness to noise and to improve the unmixing quality. We show that our method can disentangle multiple spatially overlapping labels and recover continuous maps of quantities of interest on controlled phantoms and mice experiments.