  In an effort to effectively model observed patterns in the spatial
configuration of individuals of multiple species in nature, we introduce the
saturated pairwise interaction Gibbs point process.
  Its main strength lies in its ability to model both attraction and repulsion
within and between species, over different scales.
  As such, it is particularly well-suited to the study of associations in
complex ecosystems.
  Based on the existing literature, we provide an easy to implement fitting
procedure as well as a technique to make inference for the model parameters.
  We also prove that under certain hypotheses the point process is locally
stable, which allows us to use the well-known `coupling from the past'
algorithm to draw samples from the model.
  Different numerical experiments show the robustness of the model.
  We study three different ecological datasets, demonstrating in each one that
our model helps disentangle competing ecological effects on species'
distribution.
