Multisite phosphorylation is known to be an important and dynamic mechanism for regulating the activity of transcription factors. Here we propose a stochastic pi-calculus modelling approach able to handle the complexity of post-translational modifications and to overcome the limitations of the ordinary differential equations based methods. The model can be applied without a priori assumptions to every (multisite) phosphorylation regulation for which some basic rates are known or can be indirectly set with experimental data. We apply it to the multisite phosphorylation of the Pho4 transcription factor that plays a crucial role in the phosphate starvation signalling in Saccharomyces cerevisiae using available in vitro experiments for the model tuning and validation. The innovative modelling of the sub-path with the stochastic pi-calculus allows quantitative analyses of the kinetic characteristics of the Pho4 phosphorylation, the different phosphorylation dynamics for each site (possibly combined) and the variation of the kinase activity as the reaction goes to completion. One of the performed predictions indicates that the Pho80-Pho85 kinase activity on the Pho4 substrate is nearly distributive and not semi-processive as previously found analysing only the phosphoform concentrations in vitro. This result is obtained because the model can consider and quantify the binding events without phosphorylations that cannot be experimentally measured. Thanks to the compositionality property of process algebras, we also developed the whole PHO pathway model that gives new suggestions and confirmations about its general behaviour. The potentialities of process calculi based in silico simulations for biological systems are highlighted and discussed.