A combination of molecular dynamics simulations and quantum chemistry calculations affords a good understanding of the reaction between [NPCl2](n), and the bifunctional nucleophile aminophenol HO-C6H4-NH2 in THF solution and in the presence of alkali carbonates. This reaction is regiospecific at room temperature favoring the attack by the NH2 group, and only at higher temperatures the attack by the OH is activated. It is also explained why the reaction of HO-C6H4-NH2 with {[NP(O2C12H8)](1-x)[NPCl2](x)}(n) and Cs2CO3 at reflux proceeds exclusively by the OH giving aryloxyphosphazenes with terminal NH2 substituents. Molecular dynamics simulations have been performed on model systems consisting on [NPCl2](n) chains surrounded by aminophenol molecules or ions obtained by deprotonation of either the OH or NH2 groups. Quantum mechanics calculations consisted in density functional theory (DFT) computation of the reaction mechanism of aminophenol with the phosphoranimine (NH2)Cl2P(= NH) as a simple model for the NPCl2 phosphazene units, using B3LYP approximation with 6-31G(d). Our results provide evidence to support that the macromolecular substitution reactions of poly(dichlorophosphazene) with aryloxide ions are diffusion controlled and governed by the atomic charges, while the reaction with amino-arenes are directed by hydrogen interactions between the NH, hydrogens and the phosphazene nitrogen atoms that initiate a lower activation energy stepwise mechanism competing with a concerted mechanism. The combined methodology employed in this work seems to be useful to address the regiochemistry of the macromolecular substitution on polyphosphazenes with heterobifunctional nucleophiles which is a central issue in the design of polymers with predetermined chemical composition and structure.