A method for the simulation of aggregation and breakup processes in colloidal particle suspensions is presented. The method combines a Monte Carlo algorithm to determine, on the basis of probabilistic considerations, the sequence of aggregation and breakup events, and a Discrete Element Method, built in the framework of Stokesian dynamics and contact mechanics, to accurately reproduce them. Liquid-solid suspensions subject to a uniform shear stress are investigated. The model is seen to be able to reproduce the typical dynamic steady state which is observed in colloidal suspensions under severe shearing, in which the effects of aggregation and breakup balance each other. The structural properties of the aggregates and the dynamics of the aggregation and breakup phenomena are characterized in detail. Both fragmentation and erosion are seen to contribute to the breakup process, which is characterized by an exponent similar to the one reported in the literature for compact clusters.