  In this paper, we consider the inverse problem of determining some
coefficients within a coupled nonlinear parabolic system, through boundary
observation of its non-negative solutions. In the physical setup, the
non-negative solutions represent certain probability densities in different
contexts. We innovate the successive linearisation method by further developing
a high-order variation scheme which can both ensure the positivity of the
solutions and effectively tackle the nonlinear inverse problem. This enables us
to establish several novel unique identifiability results for the inverse
problem in a rather general setup. For a theoretical perspective, our study
addresses an important topic in PDE analysis on how to characterise the
function spaces generated by the products of non-positive solutions of
parabolic PDEs. As a typical and practically interesting application, we apply
our general results to inverse problems for ecological population models, where
the positive solutions signify the population densities.
