We use the holonomic gradient method to evaluate the probability content of a simplex region under a multivariate normal distribution. This probability equals to the integral of the probability density function of the multivariate Gaussian distribution on the simplex region. For this purpose, we generalize the inclusion--exclusion identity which was given for polyhedra, to the faces of a polyhedron. This extended inclusion--exclusion identity enables us to calculate the derivatives of the function associated with the probability content of a polyhedron in general position. We show that these derivatives can be written as integrals of the faces of the polyhedron.