  We propose and analyse a novel surface finite element method that preserves
the invariant regions of systems of semilinear parabolic equations on closed
compact surfaces in $\mathbb{R}^3$ under discretisation. We also provide a
fully-discrete scheme by applying the implicit-explicit (IMEX) Euler method in
time. We prove the preservation of the invariant rectangles of the continuous
problem under spatial and full discretizations. For scalar equations, these
results reduce to the well-known discrete maximum principle. Furthermore, we
prove optimal error bounds for the semi- and fully-discrete methods, that is
the convergence rates are quadratic in the meshsize and linear in the timestep.
Numerical experiments are provided to support the theoretical findings. In
particular we provide examples in which, in the absence of lumping, the
numerical solution violates the invariant region leading to blow-up due to the
nature of the kinetics.
