Accurate extraction of physical and biochemical parameters from optoacoustic images is often impeded due to the use of unrigorous inversion schemes, incomplete tomographic detection coverage or other experimental factors that cannot be readily accounted for during the image acquisition and reconstruction process. For instance, inaccurate assumptions in the physical forward model may lead to negative optical absorption values in the reconstructed images. Any artifacts present in the single wavelength optoacoustic images can be significantly aggravated when performing a two-step reconstruction consisting in acoustic inversion and spectral unmixing aimed at rendering the distributions of spectrally-distinct absorbers. We investigate a number of algorithmic strategies with non-negativity constraints imposed at the different phases of the reconstruction process. Performance is evaluated in cross-sectional multispectral optoacoustic tomography (MSOT) recordings from tissue-mimicking phantoms and in vivo mice embedded with varying concentrations of contrast agents. Additional in vivo validation is subsequently performed with molecular imaging data involving subcutaneous tumors labeled with genetically-expressed iRFP proteins and organ perfusion by optical contrast agents. It is shown that constrained reconstruction is essential for reducing the critical image artifacts associated with inaccurate modeling assumptions. Furthermore, imposing the non-negativity constraint directly on the unmixed distribution of the probe of interest was found to maintain the most robust and accurate reconstruction performance in all experiments.