This work is the numerical analysis and computational companion of the paper by Kamgang and Sallet (Math. Biosc. 213 (2008), pp. 1-12) where threshold conditions for epidemiological models and the global stability of the disease-free equilibrium (DFE) are studied. We establish a discrete counterpart of the main continuous result that guarantees the global asymptotic stability (GAS) of the DFE for general epidemiological models. Then, we design nonstandard finite difference (NSFD) schemes in which the Metzler matrix structure of the continuous model is carefully incorporated and both Mickens' rules (World Scientific, Singapore, 1994) on the denominator of the discrete derivative and the nonlocal approximation of nonlinear terms are used in an innovative way. As a result of these strategies, our NSFD schemes are proved to be dynamically consistent with the continuous model, i.e., they replicate their basic features, including the GAS of the DFE, the linear stability of the endemic equilibrium (EE), the positivity of the solutions, the dissipativity of the system, and its inherent conservation law. The general analysis is made detailed for the MSEIR model for which the NSFD theta method is implemented, with emphasis on the computational aspects such as its convergence, or local truncation error. Numerical simulations that illustrate the theory are provided.