Axonal beading or formation of multiple beads along an axon is characteristic of many brain pathological states like Alzheimer's, Parkinson's and traumatic injuries. Despite the many existing experimental studies, the underlying mechanisms of this shape instability remain still poorly understood. In this paper, we establish a combined theoretical and numerical framework to study the governing key factors of this morphological transformation. We develop a three-dimensional (3D) non-equilibrium large deformation thermodynamic model with two main parts: the central axoplasm which is considered as a polyelectrolyte hydrogel and the encapsulating cortical membrane which is modeled as an incompressible hyperelastic layer with surface energy and growing surface. The model constitutive and evolution equations are then extracted employing thermodynamic balance principles for both bulk and surface material points. It is shown that the second law of thermodynamics indicates that the axolemma growth rate is proportional to the membrane tension which is in perfect agreement with the available experimental findings. While the developed model is general and can be extended to cover other types of axonal beadings, for the sake of simplicity, here, we focus on osmotically driven axisymmetric beadings which are compressible viscoelastic periodic modulations. We solve the corresponding governing equations using the linear perturbation method. This perturbation analysis proves that: 1) the beading instability is a rate dependent phenomenon that is controlled by the axolemma growth, 2) the initially dominant beading waves (the fastest waves) might be replaced only by longer waves which are more stable and 3) the wavelength of the fastest beads should vary roughly linearly with the axonal radius. These main findings are all in good agreement with the existing experimental results. Finally, the finite element implementation of the model is also presented to verify the results of the linear stability analysis for slow waves. The obtained axisymmetric finite element results are in good agreement with the corresponding theoretical findings.