We propose an inverse problem formulation and solution algorithm for the calibration of a multi-species biophysical model of glioblastoma (GBM) growth. The model couples GBM growth with GBM-induced deformation of the healthy parenchyma, the so called “mass effect”. The model is a multi-species partial differential equation (PDE) that models interactions between proliferative, infiltrative, and necrotic tumor cells as well as hypoxia. This PDE has several unknown parameters and fields: the pretumor brain anatomy, initial conditions, and ten coefficients representing mechanisms like diffusive proliferation, migration, and growth. A key challenge is estimating these parameters using just a single multiparametric magnetic resonance imaging (mpMRI) scan. To solve this inverse problem, we first segmented the mpMRI and then we use a single-species PDE model to invert for the tumor initial condition and the pretumor anatomy. Then we solve a second inverse problem with the multi-species PDE but without mass effect, to obtain an initial estimate of the ten unknown scalar coefficients. We use these estimates to appropriately scale the mass-effect estimates from the first inverse problem. Finally, we solve a third inverse problem with the multi-species PDE including mass effect to finalize the ten coefficient values. We solve all three inverse problems using quasi-Newton solvers, using an adjoint-based formulation for the single-species problem, and a sensitivity-based formulation for the two multi-species problems. We present preliminary results from evaluating our methodology on five subjects from the BraTS20 dataset. We show that despite the model complexity and the sparse data it is still possible to solve the inverse problem. When comparing to the single-species reconstruction, we find that the new model results in similar or better reconstruction for the overall tumor region while extracting a richer set of biophysical biomarkers.