Sorry for my late reply, we actually have a method implemented in brainstorm for statistical analysis of tPAC results based surrogate data (as Sylvain mentioned). It works with single band data (to more efficient in terms of computational cost, and interpretability of results as explained bellow). You can z-score your tPAC values with respect to surrogate data.
I would suggest not using maxPAC measure! (Could be very noisy and not really meaningful if you have a wide range for fA). Instead compute the PAC for the whole range of interest, estimate comodulogram, and then based on that find the main mode of coupling and run tPAC in the second round in a narrower band with strongest coupling using "single band" option for fA.
(As recommended in the tutorial, it would be better to select fP band also based on the your background hypothesis and power spectrum of the signal.)
It is more reasonable to run statistics on single band results (2D) both in terms of computational cost and interpretability of the result!