  We study the late time cosmological evolution of $f(R)$ theories of modified
gravity, with the matter content of the universe being that of collisional self
interacting matter. We assume that the universe is described by a flat
Friedmann-Lemaitre-Robertson-Walker metric and that it is matter and dark
energy dominated. The results of our numerical analysis for a collisional
matter $f(R)$ theory are compared with those resulting from pressure-less
matter $f(R)$ theory and from the $\mathrm{\Lambda}\mathrm{CDM}$ model. As we
shall demonstrate, the resulting picture can vary from model to model,
indicating that the effect of collisional matter in $f(R)$ theories is strongly
model dependent. Particularly, in a few cases, may give better fit to the
$\mathrm{\Lambda}\mathrm{CDM}$ model. In all studied cases, the effective
equation of state parameter does not cross the phantom divide, both in the
collisional matter and pressure-less matter $f(R)$ theories. Finally, we
thoroughly study the effects of collisional matter on one of the $f(R)$ models
that is known to provide a unified description of early time inflation and late
time acceleration. The overall picture of the evolution of the universe is not
drastically affected, apart from the matter era which is further enhanced with
an additional matter energy density term, which is of leading order. However, a
fully consistent description of the universe's evolution requires the
introduction of a dark energy compensate in the total energy density, a concept
very well known from the literature.
