We report on microscopic tight-binding modeling of surfacestates in Bi$_2$Se$_3$ three-dimensional topological insulator, based on a\textit{sp}$^3$ Slater-Koster Hamiltonian, with parameters calculated fromdensity functional theory. The effect of spin-orbit interaction on theelectronic structure of the bulk and of a slab with finite thickness isinvestigated. In particular, a phenomenological criterion of band inversion isformulated for both bulk and slab, based on the calculated atomic- andorbital-projections of the wavefunctions, associated with valence and conductionband extrema at the center of the Brillouin zone. We carry out athorough analysis of the calculated bandstructures of slabs with varyingthickness, where surface states are identified using a quantitative criterionaccording to their spatial distribution. The thickness-dependent energy gap,attributed to inter-surface interaction, and the emergence of gapless surfacestates for slabs above a critical thickness are investigated. We map out thetransition to the infinite-thickness limit by calculating explicitly themodifications in the spatial distribution and spin-character of the surfacestates wavefunction with increasing the slab thickness. Our numerical analysisshows that the system must be approximately forty quintuple-layers thick toexhibit completely decoupled surface states, localized on the oppositesurfaces. These results have implications on the effect of external perturbationson the surface states near the Dirac point.