In this paper we analyze the convergence of some commonly used Krylov subspace methods for computing the action of matrix Mittag-Leffler functions. As is well known, such functions find application in the solution of fractional differential equations. We illustrate the theoretical results by some numerical experiments.