A two-dimensional (2D) aggregation model was introduced to investigate the nucleation stage of protein crystallization. The constituent protein molecules were positioned at the lattice points of a 2D square. Four specific interaction sites were assigned to the fringe of each molecule, and through these sites the molecules are bound to each other with an effective interaction. Protein molecules were allowed to translate in the vicinity of the lattice points and to rotate within a limited angle. The partition function was calculated for aggregates consisting of 2-16 molecules. Exposed bonds, which connected exposed molecules to the rest part of the aggregates, were the weakest. The peripheral bonds, which connected molecules at the periphery of the aggregates, were weaker than the buried bonds, which were located inside the aggregates. Increasing the size of the aggregate caused both the buried and peripheral bonds to become stronger due to correlation between fluctuating bonds. The more compact an aggregate was, the more stable it was. The free energy was determined by the size, the compactness of the aggregate, and the correlation between bonds. The free-energy gain that results by generating aggregates from individual free molecules was calculated, yielding a critical size, which discriminated growing aggregates from dissociating aggregates.