Dieterich <1994> modeled the response to a stress step of a population of faults governed by rate- and state-dependent friction. This model assumes that aftershocks nucleate over areas on the fault that at the time of the main shock are already accelerating toward failure and disregards the effect of interactions among aftershocks. The main objective of this study is to examine consequences of relaxing these underlying assumptions. Aftershock activity is simulated using an inherently discrete earthquake fault model, with a fault surface governed by an approximate constitutive friction law similar to the one used by Dieterich. We find that the governing equations in nondimensional form are a function of three main parameters and explore the effect of these parameters on the simulated catalogs. We derive a simple expression for the time-dependent seismicity response to a stress step that approximates the effect of multiple interactions among aftershocks as a time-dependent stressing rate. Close match is found between the simulated seismicity response to a stress step and that predicted analytically. However, the numerical simulations show that the effect of the main shock is not only to raise the local seismicity rate but also to systematically modify the earthquake size distribution. As a result, the actual seismicity rate change early during the aftershock sequence may be higher than that predicted, whereas seismicity rate late in the sequence may be lower than that predicted. Such a modification of the earthquake size distribution can explain observations of lower b values immediately following a stress step.