Applications of fractional calculus in magnetic resonance imaging (MRI) have increased over the last twenty years. From the mathematical, computational, and biophysical perspectives, fractional calculus provides new tools for describing the complexity of biological tissues (cells, organelles, membranes and macromolecules). Specifically, fractional order models capture molecular dynamics (transport, rotation, and vibration) by incorporating power law convolution kernels into the time and space derivatives appearing in the equations that govern nuclear magnetic resonance (NMR) phenomena. Hence, it is natural to expect fractional calculus models of relaxation and diffusion to be applied to problems in NMR and MRI. Early studies considered the fractal dimensions of multi-scale materials in the non-linear growth of the mean squared displacement, assumed power-law decays of the spectral density, and suggested stretched exponential signal relaxation to describe non-Gaussian behavior. Subsequently, fractional order generalization of the Bloch, and Bloch-Torrey equations were developed to characterize NMR (and MRI) relaxation and diffusion. However, even for simple geometries, analytical solutions of fractional order equations in time and space are difficult to obtain, and predictions of the corresponding changes in image contrast are not always possible. Currently, a multifaceted approach using coarse graining, simulation, and accelerated computation is being developed to identify 'imaging' biomarkers of disease. This review surveys the principal fractional order models used to describe NMR and MRI phenomena, identifies connections and limitations, and finally points to future applications of the approach.