We report a numerical study of β-lactoglobulin aggregation using grand canonical Monte Carlo simulations of a simple lattice model in which the proteins are represented by a single lattice point and interact via a sum of a short-ranged attraction, and a long-ranged screened electrostatic repulsion. For certain values of the potential parameters we observe the so-called cluster phase, in which protein aggregates of finite size repel each other. The properties of the cluster phase are dependent on the salt concentration, the charge of the protein, and the strength of the short-ranged attraction. Disulfide bridges are modeled by covalent bonds between the lattice points, and can exchange with free thiols. Allowing the thiol–disulfide exchange leads to a severe lowering in the chemical potential of the cluster transition, or equivalently, a lower monomer density. Moreover, we find that the disulfide bridges (or the free thiol groups) are not uniformly distributed over the aggregate. The free thiol groups are significantly more abundant on the surface than in the core of the aggregate, making the surface more reactive than the inner core. This finding might explain why films made of β-lactoglobulin by cold gelation, after resolvation, reconstitute finite aggregates rather than a monomer solution.