Comparative evaluation of botulinum neurotoxin type A (BoNT/A) formulations is limited by short clinical trials and lack of long-term safety data. This study applied a validated multiscale computational model (AesthetiSIM™) to simulate 20 years of repeated treatment with eight BoNT/A formulations, integrating receptor binding kinetics, tissue diffusion, pharmacodynamic response, and HLA-weighted immunogenicity. Model parameters were calibrated against published biochemical, imaging, and randomized clinical trial data, and uncertainty was quantified through 10,000 Monte Carlo iterations. Each virtual patient underwent standardized five-site glabellar injections with clinically equivalent doses. LetibotulinumtoxinA and PrabotulinumtoxinA predicted the most balanced profiles, combining high efficacy (mean peak strain reduction 72 % and 68 %, respectively), rapid onset (2.7-3.1 days), and minimal predicted neutralizing antibody incidence over 20 years (0.4-0.6 %). DaxibotulinumtoxinA showed the longest modelled duration of effect (median 142 days), whereas AbobotulinumtoxinA and Relatox exhibited broader diffusion and higher immunogenicity risk. Overall ranking stability remained high (ρ = 0.92) across all uncertainty scenarios. These simulations align closely with available empirical data and suggest that formulation-specific differences in protein structure, diffusion behaviour, and antigenicity may contribute to observed clinical variability. Model-informed comparative analysis can complement real-world and prospective studies to optimize product selection and dosing strategies in aesthetic and therapeutic applications.