A domain decomposition method is developed for the numerical solution of nonlinear parabolic partial differential equations in any space dimension, based on the probabilistic representation of solutions as an average of suitable multiplicative functionals. Such a direct probabilistic representation requires generating a number of random trees, whose role is that of the realizations of stochastic processes used in the linear problems. First, only few values of the sought solution inside the space-time domain are computed (by a Monte Carlo method on the trees). An interpolation is then carried out, in order to approximate interfacial values of the solution inside the domain. Thus, a fully decoupled set of sub-problems is obtained. The algorithm is suited to massively parallel implementation, enjoying arbitrary scalability and fault tolerance properties. Pruning the trees is shown to increase appreciably the efficiency of the algorithm. Numerical examples conducted in 2D, including some for the KPP equation, are given. © 2009 Elsevier Inc. All rights reserved.