In recent years, a large body of the literature has been devoted to study reactive transport of solutes in porous media based on pure Lagrangian formulations. Such approaches have also been extended to accommodate second‐order bimolecular reactions, in which the reaction rate is proportional to the concentrations of the reactants. Rather, in some cases, chemical reactions involving two reactants follow more complicated rate laws. Some examples are (1) reaction rate laws written in terms of powers of concentrations, (2) redox reactions incorporating a limiting term (e.g., Michaelis‐Menten), or (3) any reaction where the activity coefficients vary with the concentration of the reactants, just to name a few. We provide a methodology to account for complex kinetic bimolecular reactions in a fully Lagrangian framework where each particle represents a fraction of the total mass of a specific solute. The method, built as an extension to the second‐order case, is based on the concept of optimal Kernel Density Estimator, which allows the concentrations to be written in terms of particle locations, hence transferring the concept of reaction rate to that of particle location distribution. By doing so, we can update the probability of particles reacting without the need to fully reconstruct the concentration maps. The performance and convergence of the method is tested for several illustrative examples that simulate the Advection‐Dispersion‐Reaction Equation in a 1‐D homogeneous column. Finally, a 2‐D application example is presented evaluating the need of fully describing non‐bilinear chemical kinetics in a randomly heterogeneous porous medium.