We study Mittag–Leffler (ML) fractional integrals involved in the solution processes of a system of coupled fractional stochastic differential equations. We introduce the ML fractional stochastic process as a ML fractional stochastic integral with respect to a standard Brownian motion. We provide some representation formulas of solution processes in terms of Mittag–Leffler fractional integrals and processes. Computable expressions of the mean functions and of the covariances of such processes are specifically given. The application in neuronal modeling is provided, and all involved functions and processes are specifically determined. Numerical evaluations are carried out and some results are shown and discussed.