We present a kinetic Monte Carlo method for simulating chemical transformations specified by reaction rules, which can be viewed as generators of chemical reactions, or equivalently, definitions of reaction classes. A rule identifies the molecular components involved in a transformation, how these components change, conditions that affect whether a transformation occurs, and a rate law. The computational cost of the method, unlike conventional simulation approaches, is independent of the number of possible reactions, which need not be specified in advance or explicitly generated in a simulation. To demonstrate the method, we apply it to study the kinetics of multivalent ligandreceptor interactions. We expect the method will be useful for studying cellular signaling systems and other physical systems involving aggregation phenomena.Proteins in cellular regulatory systems, because of their multicomponent composition, can interact in a combinatorial number of ways to generate myriad protein complexes, which are highly dynamic [1]. This feature of protein-protein interactions has been called combinatorial complexity, and it is recognized as a major barrier to understanding cell biology [1,2,3,4]. The problem of combinatorial complexity is alleviated by using a rule-based approach to model protein-protein interactions [5]. In this approach, proteins and protein complexes are represented as structured objects (graphs) and protein-protein interactions are represented as (graph-rewriting) rules that operate on these objects to modify their properties, consistent with transformations mediated by the interactions being represented. Rules can serve as definitions of individual reactions or entire reaction classes, and they can be used as generators of reactions [6,7]. The assumption underlying this modeling approach, which is consistent with the modularity of regulatory proteins [8], is that interactions are governed, at least to a first approximation, by local context that can be captured in simple rules (e.g., by the availability of binding sites on two binding partners). Rules can, in principle, be used to generate reaction networks that account comprehensively for the consequences of specified protein-protein interactions. However, the size of a rule-derived network can severely challenge conventional methods for simulating reaction kinetics [5]. For example, the rule set formulated by Danos et al. [9] implies more than 10 23 chemical species and an even greater number of reactions. It is impractical to simulate the kinetics of such a rule-derived network with the methods that are most commonly used in modeling studies of cellular regulatory systems, such as Gillespie's method [10,11]. These methods tend to be ones that are applicable in the well-mixed limit, and they are generally population based, meaning that they explicitly track populations of chemical species. The computational cost of simulation is O(log 2 M) per reaction event for efficient kinetic Monte Carlo (KMC) implementations [12,13], where M is...