A Multiscale Simulation Approach to Compute Protein–Ligand Association Rate Constants by Combining Brownian Dynamics and Molecular Dynamics

Drug-protein binding kinetic parameters are key indicators of drug efficacy, but their experimental determination is often costly and time-consuming. Computational approaches require optimization of the <i>in silico</i> approximations to achieve sufficient accuracy while remaining computationally feasible. Here, the combination of Brownian dynamics (BD) and molecular dynamics (MD) has been investigated for this purpose. BD is used for simulating long-range diffusion and diffusional encounter complex formation, whereas MD captures the subsequent formation of the bound complex, providing a detailed treatment of short-range interactions and molecular flexibility. While existing methods that employ this approach have successfully yielded estimated association rate constants (<i>k</i>_on), they often require extensive computational resources. In this work, we developed a multiscale pipeline that achieves improved computational efficiency by optimizing the sampling by BD simulation to generate an ensemble of diffusional encounter complexes in which the ligand comes very close to its protein binding site and then uses these as starting structures for MD simulation. Due to the much lower computational cost of BD simulation and the reduced MD simulation time, the approach is computationally efficient while preserving accuracy. The pipeline has been validated for a diverse set of protein-ligand complexes, varying in size, flexibility, and binding properties, yielding <i>k</i>_on values that align well with experimental measurements, as well as insights into the physical determinants of the association rate.