Voltage-gated sodium (NaV) channels generate and propagate action potentials in excitable cells, and several NaV subtypes have become important targets for pain management. The μ-conotoxins inhibit subtypes of the NaV with varied specificity but often lack of specificity to interested subtypes. Engineering the selectivity of the μ-conotoxins presents considerable complexity and challenge, as it involves the optimization of their binding affinities to multiple highly conserved NaV subtypes. In this study, a model of NaV1.4 bound with μ-conotoxin PIIIA complex was constructed using homology modeling, docking, molecular dynamic simulations and binding energy calculations. The accuracy of this model was confirmed based on the experimental mutagenesis data. The complex models of PIIIA bound with varied subtypes of NaV1.x (x = 1, 2, 3, 5, 6, 7, 8, or 9) were built using NaV1.4/PIIIA complex as a template, and refined using molecular dynamic simulations. The binding affinities of PIIIA to varied subtypes of NaV1.x (x = 1 to 9) were calculated using the Molecular Mechanics Generalized Born/Surface Area (MMGB/SA) and umbrella sampling, and were compared with the experimental values. The binding affinities calculated using MMGB/SA and umbrella sampling are correlated with the experimental values, with the former and the latter giving correlation coefficient of 0.41 (R2) and 0.68 (R2), respectively. Binding energy decomposition suggests that conserved and nonconserved residues among varied NaV subtypes have a synergistic effect on the selectivity of PIIIA.