  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.
