Standard molecular simulation methods based on classical force fields typically assume only a fixed protonation state of systems. This assumption generally only permits a limited treatment of pH effects, for example, consideration of extreme acidic/basic conditions or situations where a small number of protonation states can be explicitly enumerated. Importantly, the standard approach cannot be scaled to chemical systems with a large number of titrateable sites such as lipid head groups or assemblies of protein subunits. This talk will describe the development of a scalable and extensible method for including pH effects in molecular dynamics simulations. In contrast to other constant pH methods, the new method, based on a hybrid of non-equilibrium molecular dynamics and Monte Carlo, can be easily scaled to handle large heterogeneous systems, i.e., not only globular proteins. The applicability of the method will be discussed and motivated by previous studies in computational enzymology as well as ongoing work with membranes and membrane proteins.