Rejection-free kinetic Monte Carlo simulation of multivalent biomolecular interactions

The system-level dynamics of multivalent biomolecular interactions can be simulated using a rule-based kinetic Monte Carlo method in which a rejection sampling strategy is used to generate reaction events. This method becomes inefficient when simulating aggregation processes with large biomolecular complexes. Here, we present a rejection-free method for determining the kinetics of multivalent biomolecular interactions, and we apply the method to simulate simple models for ligand-receptor interactions. Simulation results show that performance of the rejection-free method is equal to or better than that of the rejection method over wide parameter ranges, and the rejection-free method is more efficient for simulating systems in which aggregation is extensive. The rejection-free method reported here should be useful for simulating a variety of systems in which multisite molecular interactions yield large molecular aggregates.

💡 Research Summary

The paper addresses a fundamental bottleneck in the simulation of multivalent biomolecular interactions using rule‑based kinetic Monte Carlo (KMC) methods. Traditional implementations rely on a rejection‑sampling scheme: a reaction rule is chosen at random, a candidate pair of binding sites is drawn, and the event is accepted only if the pair is actually compatible. While mathematically exact, this approach becomes increasingly inefficient when the system forms large aggregates, because the number of possible site‑pair combinations grows dramatically and most random draws are rejected.

To overcome this limitation, the authors propose a rejection‑free algorithm that directly samples reaction events according to their exact probabilities without any “trial‑and‑error” step. The key technical innovations are: (1) a dynamic partner‑list for each molecule and each binding site that records all currently available interaction partners; (2) a data structure combining hash tables and linked lists that allows constant‑time insertion, deletion, and lookup of partners as bonds form or break; and (3) a weighted‑selection step in which the total propensity ΣR of all possible reactions is computed, each individual reaction i is assigned a weight Ri/ΣR, and a single uniform random number determines which reaction fires. This is essentially the Gillespie direct method applied at the level of individual binding events, guaranteeing statistical exactness while eliminating the costly rejection loop.

The authors provide a theoretical complexity analysis. In the rejection‑free scheme, selecting an event and updating the partner lists requires O(1) to O(log N) operations, where N is the size of the largest aggregate. By contrast, the rejection‑based method can degrade to O(N) or worse because each trial may need to scan a growing list of potential partners. Consequently, the new algorithm’s runtime is essentially independent of aggregate size, making it well suited for highly polymerizing systems.

Two benchmark models are used for validation. Model 1 consists of monovalent ligands interacting with multivalent receptors; Model 2 features two multivalent species that can cross‑link. For each model the authors sweep binding (k_on) and unbinding (k_off) rates, initial concentrations, and the number of binding sites per molecule. Across a broad parameter space, the rejection‑free method matches the rejection‑based method in terms of the statistical distributions of bound complexes, mean aggregate size, and time‑course kinetics, confirming that no bias is introduced. Performance measurements show that the new method is typically 20 %–70 % faster for moderate aggregation and up to five‑fold faster when extensive aggregation occurs (high concentration, low k_off).

The paper’s contributions can be summarized as follows: (i) a rigorously exact, rejection‑free KMC framework tailored to multivalent interactions; (ii) an implementation strategy that uses lightweight data structures to maintain up‑to‑date partner lists with minimal overhead; (iii) extensive benchmarking that demonstrates superior scalability in regimes where traditional methods become prohibitive. The authors also discuss limitations: the current implementation assumes a fixed set of binding sites per molecule, and extending it to systems where sites appear or disappear dynamically (e.g., conformational changes) would require additional bookkeeping. Memory consumption and parallelization are identified as future engineering challenges for simulations involving millions of molecules.

Finally, the authors outline several promising directions: integrating the rejection‑free engine with GPU‑accelerated parallel KMC, coupling it to spatial diffusion models, and applying it to biologically realistic scenarios such as immune‑cell receptor clustering, phase‑separated signaling condensates, and synthetic biomaterials where multivalent cross‑linking drives macroscopic properties. By removing the rejection bottleneck, this work paves the way for efficient, high‑fidelity simulations of complex multivalent networks that were previously out of reach.