If you take an infinite 2D system (so no boundaries) and find the eigenstates of the system (Bloch waves) you can determine, say, the dispersion relation as a function of the wave-vector, \epsilon_{\mathbf{k}}, which is a continuous variable in the Brillouin zone. If you now do the same analysis for a finite sized system, you'll see that \epsilon_{\mathbf{k}} assumes exactly the same form (as determined by the lattice), except that \mathbf{k} can now assume only discrete values within the Brillouin zone. As you tend the size of the system to infinity (when the length is much larger than the lattice constant) the discrete spectrum starts to look continuous (the spacing between allowed values of \mathbf{k} becomes very small). The same applies for periodic boundary conditions. As f95toli pointed out, the boundaries play a lesser role for large systems, but become important for "finite" sized systems.
I suppose mathematically, it would be related to the fact that you can turn the following sum into an integral when the system size, L^2 becomes very large:
<br />
\frac{1}{L^2}\sum_{\mathbf{k}\in BZ} \to \frac{1}{\Omega}\int_{\Omega} d^2\mathbf{k}<br />
(\Omega is the area of the BZ. You usually end up evaluating such summations - in the calculation of the Green's function for example, and other quantities).
