A multiscale simulation approach to compute protein-ligand association rate constants by combining Brownian dynamics and molecular dynamics

Drug-protein binding kinetic parameters are a key indicator of drug efficacy in vivo, but their experimental determination is often costly and time-consuming. Consequently, computational approaches are becoming an alternative to experimental assays. However, in silico approximations need optimization to achieve sufficient accuracy while remaining computationally feasible. Multiscale methods combining Brownian dynamics (BD) and molecular dynamics (MD) stand out as an effective solution. BD is used for simulating long-range diffusion and diffusional encounter complex formation, while MD captures the subsequent formation of the bound complex with a detailed treatment of short-range interactions and molecular flexibility. While existing methods that employ this multiscale approach have successfully yielded estimated association rate constants (kon), they often require extensive computational resources. In this work, we developed a multiscale workflow that improves efficiency over previously reported methods by optimizing the sampling by BD simulation to generate and ensemble of diffusional encounter complexes in which the ligand comes very close to the active site and then using these as starting structures for MD simulation. Due to the much lower computational costs of BD simulation and the reduced MD simulation time, the approach is computationally efficient while preserving accuracy. The pipeline is validated on a diverse set of protein-ligand complexes, varying in size, flexibility, and binding properties, yielding kon values that align well with experimental data.