The mean field limits of systems of interacting diffusions (also called stochastic interacting particle systems (SIPS)) have been intensively studied since McKean \cite{mckean1966class}. The interacting diffusions pave a way to probabilistic representations for many important nonlinear/nonlocal PDEs, but provide a great challenge for Monte Carlo simulations. This is due to the nonlinear dependence of the bias on the statistical error arising through the approximation of the law of the process. This and the fact that particles/diffusions are not independent render classical variance reduction techniques not directly applicable and consequently make simulations of interacting diffusions prohibitive. In this article, we provide an alternative iterative particle representation, inspired by the fixed point argument by Sznitman \cite{sznitman1991topics}. This new representation has the same mean field limit as the classical SIPS. However, unlike classical SIPS, it also allows decomposing the statistical error and the approximation bias. We develop a general framework to study integrability and regularity properties of the iterated particle system. Moreover, we establish its weak convergence to the McKean-Vlasov SDEs (MVSDEs). One of the immediate advantages of iterative particle system is that it can be combined with the Multilevel Monte Carlo (MLMC) approach for the simulation of MVSDEs. We proved that the MLMC approach reduces the computational complexity of calculating expectations by an order of magnitude. Another perspective on this work is that we analyse the error of nested Multilevel Monte Carlo estimators, which is of independent interest. Furthermore, we work with state dependent functionals, unlike scalar outputs which are common in literature on MLMC. The error analysis is carried out in uniform, and what seems to be new, weighted norms.