Computer-implemented methods, systems, and computer-readable

Abstract

The invention relates to methods and systems for computing an alchemical free-energy difference, the method comprising : a) obtaining at least one collective variable describing the molecular system b) during the simulation, updating on the fly and applying to the alchemical parameter λ a biasing scheme comprising: an adaptive biasing force configured to drive sampling of collective variable toward a flattened distribution, and a bias potential added as an additional potential energy term, the bias potential being defined as a function of a target probability distribution of collective variable, wherein the bias potential is constructed using Gaussian kernels; and c) computing the alchemical free-energy difference by thermodynamic integration of energy derivatives calculated from the collected information with respect to the collective variable obtained during the step b).

Description

[0001] COMPUTER-IMPLEMENTED METHODS, SYSTEMS, AND COMPUTER-READABLE MEDIA FOR COMPUTING FREE-ENERGY SURFACES IN MOLECULAR SYSTEMS

[0002] Field of the invention

[0003] The present disclosure relates to computer-implemented molecular simulation methods and systems for computing free-energy surfaces of molecular systems. The invention has several applications such as in drug discovery, molecular dynamics, and protein folding studies.

[0004] Description of related art

[0005] Molecular simulation techniques, and in particular molecular dynamics (MD), are used to model the structure, dynamics, and thermodynamics of molecular systems, including solvated molecules, biomacromolecules, and molecular complexes.

[0006] One of the key challenges in these simulations is the accurate computation of free-energies, which is essential for predicting the binding affinities of ligands to their targets.

[0007] Free-energy surfaces (and free-energy differences derived therefrom) are used in computational chemistry for tasks including, for example, the estimation of binding affinities and absolute binding free energies. In practice, molecular simulations may explore relevant collective variables slowly due to free-energy barriers, local minima, and limited ergodicity. Classical methods for computing free-energies, such as Free Energy Perturbation (FEP) and Thermodynamic Integration (Tl), are widely used but often suffer from slow convergence and high computational costs. These methods require extensive sampling to achieve accurate results, which can be time-consuming and resource-intensive.

[0008] Enhanced sampling techniques, such as Adaptive Biasing Force (ABF) and Metadynamics, have been developed to accelerate the convergence of molecular dynamics simulations. These techniques aim to overcome free-energy barriers and improve the sampling of rare events, leading to more efficient computations of free-energies.

[0009] To address these challenges, several enhanced sampling algorithms have been developed. Notably, Collective Variable (CV)-based importance-sampling techniques such as Umbrella Sampling (US) (Torrie, G.M. et al., 1977), Adaptive Biasing Force (ABF) (Darve, E. et al., 2001 ; Darve, E. et al., 2008; Comer, J. et al., 2015), MetaDynamics (MtD) (Laio, A. et al., 2002; Barducci, A. et al., 2011), its recent evolution, On-the-fly Probability Enhanced Sampling (OPES) (Invernizzi, M. et al., 2020; Invernizzi M. et al., 2022), and Temperature-Accelerated Molecular Dynamics (TAMD) (Maragliano , L. et al., 2006) have demonstrated significant success. These algorithms rely on the definition of CVs, reduced dimensions along which biasing forces, potentials, or probabilities are applied to facilitate sampling of the system’s configuration space. Each enhanced-sampling method has its advantages andlimitations (Henin, J. et al., 2022). To overcome these challenges, hybrid methods have been developed to combine different techniques and mitigate their individual limitations.

[0010] Another class of methods, called alchemical approaches, estimates free energy differences associated to unphysical “alchemical” changes of the molecular Hamiltonian by scaling some key interactions with an alchemical parameter Ae [0; 1]. They have proven useful to compute solvation as well as binding free energies. Traditional techniques involve several simulations at fixed values and the reconstruction of the free energy difference of interest through an estimator such as the Free Energy Perturbation (FEP) or the Thermodynamic Integration (Tl) one.

[0011] Alternatively, another approach called lambda-dynamics was introduced, in which the coupling parameter is treated as a dynamical variable through an extended Lagrangian / Hamiltonian scheme. Building upon this concept, we introduced a new alchemical method, Lambda-Adaptive Biasing Force (Lambda-ABF) (Lagardere et al., 2024), that leverages lambda dynamics in combination with multiple-walker ABF, enabling efficient sampling of as a collective variable (CV). It has been shown to be robust and to improve sampling efficiency compared to standard fixed-lambda methods(Lagardere et al. 2024). Additionally, in the context of binding simulations, it leverages Distance-to-Bound-Configuration (DBC) (Salari et al., 2018) restraints to keep the ligand within the binding pocket, limiting its translational, rotational, and conformational fluctuations. Compared to fixed-lambda methods, Lambda-ABF not only reduces the computational cost but is also arguably simpler as it bypasses the definition of a schedule and portable thanks to its implementation within the CoIvars library (Fiorin et al. 2013, Fiorin et al. 2024). ABF applies biasing forces along the transition coordinate to flatten the sampled free energy

[0012] landscape toward a uniform distribution. It is well understood mathematically (Comer et al., 2015) and is associated to a local (unconstrained) Tl free-energy estimator. A key limitation arises with barriers along other degrees of freedom which can lead to kinetic trapping in some regions and overall non-ergodic sampling and also in the diffusive regime when (local) convergence is reached.

[0013] Hence, despite these advancements, challenges remain in developing methods that can efficiently compute free-energies while maintaining accuracy and reducing computational load. Addressing these challenges could lead to faster and more effective simulations in various applications, including drug discovery and molecular dynamics studies.

[0014] Summary of the invention

[0015] The following sets forth a simplified summary of selected aspects, embodiments and examples of the present invention for the purpose of providing a basic understanding of the invention. However, this summary does not constitute an extensive overview of all theaspects, embodiments and examples of the invention. Its sole purpose is to present selected aspects, embodiments and examples of the invention in a concise form as an introduction to the more detailed description of the aspects, embodiments and examples of the invention that follow the summary. The invention aims to overcome the disadvantages of the prior art. The present invention provides a computer-implemented method for computing free-energies of molecular systems using enhanced sampling techniques.

[0016] The invention relates to a computer-implemented method for computing an free-energy difference between two end-state Hamiltonians of a molecular system by molecular dynamics simulation, the method comprising:

[0017] - a) obtaining at least one collective variable describing the molecular system; - b) During the molecular dynamics simulation, applying a biasing scheme to the at least one collective variable, said a biasing scheme being updated on the fly, said biasing scheme comprising:

[0018] o an adaptive biasing force configured to drive sampling of the least one collective variable toward a flattened distribution, and

[0019] o a bias potential added as an additional potential energy term, the bias potential being defined as a function of a target probability distribution of the least one collective variable, wherein the target probability distribution is a well-tempered distribution, and wherein the bias potential is constructed using Gaussian kernels; and

[0020] - c) computing the free-energy difference by thermodynamic integration of energy derivatives calculated from the collected information, with respect to at least one collective variable, obtained during the step b).

[0021] The invention relates to a computer-implemented method for computing free-energy surface in a molecular system, the method comprising:

[0022] - Step a) : Obtaining one or several collective variable(s) (CVs) describing the molecular system; Preferably step a) also comprise defining one or several collective variable(s) (CVs) describing the molecular system;

[0023] - Step b) : Applying an adaptive biasing force to dynamically sample the collective variable (CV) of the molecular system under a flattened distribution and collecting information about the molecular system’s energy; said biasing force being updated on the fly and including the use of a bias potential expressed as a function of a target probability distribution of the molecular system; and

[0024] - Step c) : Performing thermodynamic integration of energy derivatives calculated from the collected information during the step b) to compute the free-energy surface in the molecular system.The invention relates to a computer-implemented method for computing binding free-energy surface in a molecular system, the method comprising:

[0025] - Step a) : Obtaining one or several collective variable(s) (CVs) describing the molecular system; Preferably step a) also comprise defining one or several collective variable(s) (CVs) describing the molecular system;

[0026] - Step b) : Applying an adaptive biasing force to dynamically sample the alchemical lambda parameter (A) of the molecular system under a flattened distribution and collecting information about the molecular system’s energy; said biasing force including the use of a bias potential expressed as a function of a target probability distribution of the molecular system; and

[0027] - Step c) performing thermodynamic integration of energy derivatives calculated from the collected information during the step b) to compute the alchemical binding free-energy surface in the molecular system.

[0028] The invention relates preferably to a computer-implemented method for computing an alchemical free-energy difference between two end-state Hamiltonians of a molecular system by molecular dynamics simulation, the method comprising :

[0029] - a) obtaining at least one collective variable describing the molecular system, said at least one collective variable comprising an alchemical parameter A that interpolates between the two end-state Hamiltonians;

[0030] - b) during the molecular dynamics simulation, updating on the fly and applying a biasing scheme to the alchemical parameter A, said biasing scheme comprising:

[0031] o an adaptive biasing force configured to drive sampling of A toward a flattened distribution, and

[0032] o a bias potential added as an additional potential energy term, the bias potential being defined as a function of a target probability distribution of A, wherein the target probability distribution is a well-tempered distribution, and wherein the bias potential is constructed using Gaussian kernels; and

[0033] - c) computing the alchemical free-energy difference by thermodynamic integration of energy derivatives calculated from the collected information with respect to A obtained during the step b).

[0034] In such methods, an alchemical parameter A scales selected interactions in the molecular Hamiltonian to connect two end states and enable computation of free energies, like binding free energies. Traditional fixed-A approaches typically require multiple simulations at discrete A values and subsequent reconstruction of the free-energy difference. In contrast, lambda-dynamics treats A as a dynamical variable, which can reduce computational cost and simplify execution as described in the examples. However, biasing a single coordinate such as A can still lead to inefficiencies when barriers in other degrees of freedom cause kinetic trapping and non-ergodic sampling. To address this technical limitation, the invention provides a hybrid computer-implemented method that combines (i) an Adaptive Biasing Force scheme applied to A (Lambda-ABF) and (ii) a MetaDynamics-family bias potential based on a target probability distribution (OPES-type, well-tempered) constructed using Gaussian kernels. The hybrid approach exploits complementary effects: ABF reduces free-energy barriers along A and supplies a local Tl estimator, while the kernel-based target-distribution bias fills energy valleys via a history-dependent potential term, thereby pushing the system out of local minima and mitigating kinetic trapping. As illustrated in examples, this combination improves convergence speed by up to about nine times compared with Lambda-ABF alone, while maintaining high accuracy, and reducing computational resource requirements, enabling large-scale / high-throughput applications.

[0035] Hence, although OPES-Explore alone yields inconsistent AG results across repeats for alchemical legs, the use of OPES-Explore in a qualitatively different regime, i.e. acting mainly as a controlled “push-out of local minima” with AE not tied to the free-energy barrier, while A-ABF provides the primary driving and rapid Tl convergence, results in up to nine-fold acceleration with maintained accuracy.

[0036] The synergistic effect of this combination is not obvious as ABF / TI and OPES solve enhanced sampling in fundamentally different ways. ABF / TI relies on mean forces / derivatives and Tl for recovery of free energy. While OPES relies on target-distribution probability shaping and self-consistent probability estimation (weighted KDE, etc.).

[0037] Introducing a second, target-based bias potential changes the sampling distribution and must be done in a way that remains consistent with the computed derivatives and the intended output.

[0038] Accordingly, the invention provides computational methods and simulation protocols that (i) accelerate convergence of alchemical free-energy calculations, particularly for absolute binding free energies, (ii) reduce sensitivity to kinetic trapping and slow orthogonal relaxation, (iii) minimize user intervention (e.g., in selecting A schedules and tuning numerous hyperparameters), and (iv) remain compatible with computationally demanding but higher-accuracy interaction models, including polarizable force fields, while preserving the thermodynamic rigor required for reliable binding affinity prediction.

[0039] It is worth noting that the output is specifically adapted for binding affinity workflows, and the method yields computational resource reduction. Hence, the output (e.g. free-energy surface) is not merely presented for cognitive evaluation; it is obtained through thermodynamic integration of derivatives gathered during biased sampling and is a functionaloutput that is specifically adapted to technical use in molecular simulation workflows (e.g., to drive molecular design decisions or to validate binding energetics).

[0040] According to other optional features of the method according to the invention, it can optionally include one or more of the following characteristics alone or in combination:

[0041] - at least one collective variable is selected among alchemical variable, distance, angle, dihedral, root mean square deviation, coordination number and / or machine-learning based collective variables.

[0042] - the at least one collective variable comprise an alchemical parameter A that interpolates between the two end-state Hamiltonians and the biasing scheme is applied to the alchemical parameter A; the adaptive biasing force being configured to drive sampling of the alchemical parameter A toward a flattened distribution and the bias potential being added as an additional potential energy term, the bias potential being defined as a function of a target probability distribution of the alchemical parameter A.

[0043] - there is only one collective variable. Choosing a single CV enables convergence with bounded cost, while accuracy / convergence remains acceptable.

[0044] - there is several collective variables.

[0045] - at the step b), a weighted kernel density estimation is employed. A weighted KDE yields measurable convergence gains.

[0046] - the weighted kernel density estimation uses an adaptive kernel width that decreases as the number of effective samples increases. This can provide a coarse-grained bias potential in early simulation stages and progressively higher resolution as more data accumulates. This can improve bias accuracy and stability without requiring manual parameter tuning.

[0047] - at the step b), it comprises a step of estimating an underlying probability distribution of the CV and a step of adjusting the bias potential so that the simulated distribution of the CVs approaches the target probability distribution. This can lead to faster and more reliable convergence.

[0048] - estimating the underlying probability distribution comprises performing a weighted kernel density estimation of the CV distribution using simulation samples, weighting each sample by a factor related to the current bias to approximate the unbiased probability; preferably it also comprises compressing the kernel representation by merging kernels in previously sampled regions to limit computational growth. Hence, the distribution estimate is refined over time with bounded memory usage.- at the step b), the adaptive biasing force is computed from an on-the-fly estimate of the equilibrium (unbiased) target probability distribution of the molecular system.

[0049] - at the step b), the adaptive biasing force is computed from a biasing potential issued from an on-the-fly estimate of the current biased probability distribution of the molecular system. This uses the probability density of CV observed under the bias as the basis to update the bias rather than the target distribution (equilibrium distribution).

[0050] - at the step b), the adaptive biasing force is computed from a biasing potential issued from an on-the-fly estimate of the well-tempered distribution of the molecular system being sampled. This improves the memory usage and runtime of the method.

[0051] it is a computer-implemented method for computing alchemical binding free- energy in the molecular system.

[0052] - The alchemical parameter A is between 0 and 1. This ensures valid interpolation between end states, avoids invalid Hamiltonians and numerical instabilities and hence prevents wasted compute.

[0053] it comprises a step of defining reflective boundary conditions for A=0 and A=1 , ensuring that said alchemical parameter remains within [0,1] during the simulation.

[0054] - the Gaussian kernels are defined with parameters comprising at least one from:

[0055] Gaussian kernel width and barrier parameter; preferably barrier parameter.

[0056] - the barrier parameter (AE) is selected to allow efficient transitions between basins while preventing access to irrelevant high-energy states.

[0057] - an adaptive sigma is used to determine the Gaussian kernel width.

[0058] - a frequency for kernel deposition is set as at least one of: a number of simulation steps and a deposition time interval.

[0059] - the step b) is applied concurrently by multiple-walker sampling instances of the alchemical parameter A.

[0060] - the multiple-walker sampling instances exchange or share adaptive biasing force data so as to collectively explore different regions of phase space.

[0061] - each multiple-walker sampling instances periodically transmits locally accumulated biasing data to a central data structure or to all other multiple-walker sampling instances.

[0062] - each multiple-walker sampling instances :

[0063] updates the molecular system’s atomic coordinates and velocities, updates the alchemical lambda parameter A,- calculates a local mean force on alchemical lambda parameter A and adjusts the adaptive bias, and

[0064] periodically exchanges mean force or bias information with other multiplewalker sampling instances;

[0065] o aggregating the simulation data from all multiple-walker sampling instances to compute a free-energy difference for each value of A; and o integrating the aggregated mean force data to obtain the alchemical free-energy of the molecular system.

[0066] - each multiple-walker sampling instances :

[0067] updates the molecular system’s atomic coordinates and velocities, updates the alchemical lambda parameter A,

[0068] - calculates a local mean force on alchemical lambda parameter A and adjusts the adaptive bias including the use of a bias potential expressed as a function of a target probability distribution of the molecular system, and periodically exchanges mean force or bias information with other multiplewalker sampling instances;

[0069] o aggregating the simulation data from all multiple-walker sampling instances to compute a free-energy gradient for each value of A; and o integrating the aggregated mean force data to obtain the alchemical free-energy of the molecular system.

[0070] - the molecular system comprises a ligand and A interpolates between a fully coupled state and a decoupled state of the ligand.

[0071] - a binding free energy is calculated by separately decoupling van der Waals interactions and electrostatic interactions.

[0072] - polarizabilities and permanent multipoles of the ligand are scaled down to 0 in electrostatic legs, and van der Waals interactions between atoms of the ligand and all other atoms are scaled down to 0 in van der Waals legs.

[0073] - the van der Waals legs leverage softcore interactions.

[0074] it further comprises leveraging Distance-to-Bound-Configuration (DBC) restraints to at least a part of the molecular system; preferably a ligand is partially restrained by a distance-to-bound-configuration (DBC) variable during equilibration.

[0075] - the distance-to-bound-configuration (DBC) variable is defined as a root mean square deviation of ligand atoms for each frame, with an alignment performed relative to atoms of a receptor binding site.- a DBC cutoff is selected as a DBC value within a 95% interval of a DBC distribution monitored during plain molecular dynamics of the ligand, and the selected DBC cutoff is used as a cutoff for the DBC restraint.

[0076] - a flat-bottomed harmonic restraint is applied to the DBC above the DBC cutoff. - a free energy cost associated with a release of the DBC restraint is computed in a gas phase through thermodynamic integration by progressively releasing the DBC restraint to a compatible harmonic distance restraint which is computed analytically.

[0077] - selecting ligand atoms used for the DBC variable comprises monitoring a root mean square fluctuation of ligand heavy atoms over a plain molecular dynamics simulation and selecting atoms based on the root mean square fluctuation.

[0078] - A is a coupling parameter controlling the interpolation between two end-state Hamiltonians (e.g., a ligand fully non-interacting vs fully interacting).

[0079] - performing thermodynamic integration comprises continuously integrating an ensemble-averaged derivative of the molecular system’s potential energy with respect to the alchemical parameter A over a range of A values to obtain the alchemical free-energy of the molecular system.

[0080] - the molecular simulation employs a force field selected from: fixed-charge force fields (e.g. AMBER, CHARMM, OPLS-AA), machine learning force fields, coarsegrained force fields, reactive force fields, bond-order based force fields, ionic force fields and / or polarizable force fields (e.g. AMOEBA); or the molecular simulation employs an Hamiltonian obtained through Quantum Mechanical / Molecular Mechanical (QM / MM) methods, and the lambda-dynamics adaptive biasing force approach is applied under either one or several types of force-field models.

[0081] - the thermodynamic integration comprises numerically integrating mean forces exerted along said alchemical parameter A, thereby obtaining a free-energy difference between a bound state and an unbound state.

[0082] - free-energy is monitored in real time by evaluating partial integrals of a mean force with respect to A and analyzing resulting estimates over cumulative simulation time.

[0083] it comprises the computation of potential derivatives with respect to the alchemical lambda parameter through interpolation of polarization between the end states which gives immediately the associated derivatives as the difference of these (for the polarizable force field).

[0084] - step (c) comprises computing a free-energy profile as a function of A by integrating said mean forces with respect to A, and obtaining the alchemical free-energy difference between the two end-state Hamiltonians from said integrated profile.

[0085] - step (c) is carried out during step (b) so as to provide an on-the-fly estimate of the alchemical free-energy difference during the simulation.

[0086] - step (c) comprises monitoring free-energy in real time by evaluating partial integrals of the mean force with respect to A and analyzing resulting estimates over cumulative simulation time.

[0087] - step (c) comprises integrating mean forces over a range of A values sampled continuously in a single A-dynamics simulation, without breaking the A range into discrete fixed-A windows.

[0088] - step (c) comprises aggregating mean-force data obtained from a plurality of multiple-walker sampling instances and integrating the aggregated mean-force data with respect to A to obtain the alchemical free energy.

[0089] - aggregating comprises computing a free-energy gradient (mean-force profile) for each value of A from the plurality of multiple-walker sampling instances, and integrating said free-energy gradient with respect to A.

[0090] - a standard free energy of binding is determined using a thermodynamic cycle based on (i) a complex phase in which the ligand is decoupled in complex with a protein and (ii) a solvent phase in which the ligand is decoupled in bulk solvent, each phase’s decoupling free energy being computed by step (c).

[0091] - the mean forces integrated in step (c) comprise potential-energy derivatives with respect to A computed during the simulation.

[0092] the molecular simulation employs a polarizable force field, and wherein the potential-energy derivatives comprise a derivative contribution associated with polarization energy.

[0093] it is used to: compute a binding free energy of a ligand to a target; compute an absolute binding free energy of a ligand; compute a binding affinity between a ligand and a receptor; compute a free-energy profile for protein folding; compute free-energy differences between molecular conformational states; compute activation free energies; compute reaction free energies; compute solvation free energies; compute partition coefficients; compute dissociation constants; compute association constants; compute relative stabilities of molecular conformers; compute free-energy differences in enzyme-substrate complexes; compute free energies for conformational transitions; compute free energies for phase transitions; compute free energies of adsorption; compute free energies for ligand- induced conformational changes; compute free energies in molecular docking simulations; compute free energies for applications in drug design; compute freeenergies for stability prediction of molecular systems; compute free energies in polymer systems; compute free energies in solid-state systems; compute free energies in solvated systems; compute free energies for protein-ligand interactions; compute free energies for nucleic acid-ligand interactions; compute free energies in molecular dynamics simulations; compute free energies for conformational sampling; and / or compute free energies for estimating reaction rates.

[0094] According to another aspect, the invention relates to a system for computing free-energy surface in a molecular system, preferably for computing an alchemical free-energy difference between two end-state Hamiltonians of a molecular system by molecular dynamics, the system comprising one or several processors and one or several non-transitory computer-readable memories storing instructions which, when executed by the processor(s), cause the processor(s) to perform the method according to the invention.

[0095] In particular, the invention relates to a system for computing an alchemical free-energy difference between two end-state Hamiltonians of a molecular system by molecular dynamics, the system comprising one or several processors and one or several non-transitory computer-readable memories storing instructions which, when executed by the processor(s), cause the processor(s) to perform the following steps:

[0096] - a) obtaining at least one collective variable describing the molecular system, said at least one collective variable comprising an alchemical parameter A that interpolates between the two end-state Hamiltonians;

[0097] - b) during the simulation, updating on the fly and applying to the alchemical parameter A a biasing scheme comprising:

[0098] o an adaptive biasing force configured to drive sampling of A toward a flattened distribution, and

[0099] o a bias potential added as an additional potential energy term, the bias potential being defined as a function of a target probability distribution of A, wherein the target probability distribution is a well-tempered distribution, and wherein the bias potential is constructed using Gaussian kernels; and (c) computing the alchemical free-energy difference by thermodynamic integration of energy derivatives calculated from the collected information with respect to A obtained during the step b).

[0100] According to another aspect, the invention relates to a non-transitory computer-readable medium storing instructions which, when executed by a processor, cause the processor to perform a computer-implemented method for computing a free-energy surfaceof a molecular system, preferably for computing an alchemical free-energy difference between two end-state Hamiltonians of a molecular system by molecular dynamics, according to the invention.

[0101] In particular, the invention relates to a non-transitory computer-readable medium storing instructions which, when executed by a processor, cause the processor to perform a computer-implemented method for computing a free-energy surface of a molecular system, preferably for computing an alchemical free-energy difference between two end-state Hamiltonians of a molecular system by molecular dynamics, the method comprising:

[0102] - a) obtaining at least one collective variable describing the molecular system, said at least one collective variable comprising an alchemical parameter A that interpolates between the two end-state Hamiltonians;

[0103] - b) during the simulation, updating on the fly and applying to the alchemical parameter A a biasing scheme comprising:

[0104] o an adaptive biasing force configured to drive sampling of A toward a flattened distribution, and

[0105] o a bias potential added as an additional potential energy term, the bias potential being defined as a function of a target probability distribution of A, wherein the target probability distribution is a well-tempered distribution, and wherein the bias potential is constructed using Gaussian kernels; and (c) computing the alchemical free-energy difference by thermodynamic integration of energy derivatives calculated from the collected information with respect to A obtained during the step b).

[0106] Brief description of the drawings

[0107] The foregoing and other objects, features and advantages of the present invention will become more apparent from the following detailed description when taken in conjunction with the accompanying drawings in which:

[0108] FIG. 1 is a flowchart of an overall free-energy surface computation method using on-the-fly biasing and thermodynamic integration.

[0109] FIG. 2 illustrates an example structural representation of a BRD4-ligand complex, (a) Cartoon representation of the BRD4 bromodomain structure in complex with ligand 1 (PDB ID: 4OGJ), with the BC-loop highlighted using a dashed marker, (b) Binding mode of BRD4 with the ligand, highlighting key interacting residues Asn140 and Tyr97 in cyan.

[0110] Crystallographically observed water molecules are shown as ball-and-stick representations, along with the hydrogen bond network bridging the ligand and protein viaTyr97 and Gln85 residues. The protein is represented as a transparent cartoon, (c) The 2D structures of thecompounds analyzed with this invention are shown and labeled with Arabic numerals in descending order of binding affinity.

[0111] FIG. 3 illustrates example free-energy surfaces (FES) of a X2 dihedral angle for an Asn140 residue obtained from plain OPES simulations starting from an Apo structure. (PDB ID:

[0112] 4LYI). Two different rotamers of Asn140 corresponding to each minimum are represented. The protein is depicted in a cartoon representation in yellow, while Asn140 is shown as a stick representation. X2 is the torsion angle between the nitrogen and the carbonyl carbon of the amide group. The transparent regions indicate the associated errors of the FES, calculated over block analysis.

[0113] FIG. 4 illustrates an example comparison of experimental versus calculated AG values. In one example, calculated AG values and errors are a mean and a standard error from the mean from three repeats for each ligand. In one example, a dark shaded region spans ± 1 kcal / mol and a lighter region spans ± 2 kcal / mol. In one example, a color bar represents an absolute value difference between experimental and computed values, and Pearson's r, RMSE, and MAE are 0.81, 1.10, and 0.90 kcal / mol, respectively.

[0114] FIG. 5 illustrates example AG convergence comparison between Lambda-ABF and Lambda-ABF-OPES. In one example, AG over time is shown for an electrostatic (ELE) leg and a van der Waals (VDW) leg of a ligand, calculated using Lambda-ABF and Lambda-ABF-OPES methods, and each panel includes a zoomed-in plot of a converged area.

[0115] FIG. 6 illustrates example AG convergence comparison between Lambda-ABF-OPES and Lambda-ABF-WTMD. In one example, AG over time is shown for an electrostatic (ELE) leg and a van der Waals (VDW) leg of a ligand, calculated using Lambda-ABF-OPES and Lambda-ABF-WTMD methods, respectively, and each panel includes a zoomed-in plot of a converged area.

[0116] FIG. 7 illustrates an example FES of a VDW leg in a solvent phase of a ligand across multiple replicas.

[0117] FIG. 8 illustrates atoms Involved in the Definition of DBC for Restraint. Atoms chosen from the ligand are highlighted in orange, and atoms from the protein are shown as yellow spheres.

[0118] Several aspects of the present invention are disclosed with reference to flow diagrams and / or block diagrams of methods, devices and computer program products according to embodiments of the invention. On the figures, the flow diagrams and / or block diagrams show the architecture, the functionality and possible implementation of devices or systems or methods and computer program products, according to several embodiments of the invention. For this purpose, each box in the flow diagrams or block diagrams may represent a system, a device, a module or code which comprises several executable instructions for implementing the specified logical function(s). In some implementations, the functionsassociated with the box may appear in a different order than indicated in the drawings. For example, two boxes successively shown, may be executed substantially simultaneously, or boxes may sometimes be executed in the reverse order, depending on the functionality involved. Each box of flow diagrams or block diagrams and combinations of boxes in flow diagrams or block diagrams may be implemented by special systems that perform the specified functions or actions or perform combinations of special equipment and computer instructions.

[0119] Detailed description

[0120] Hereinafter, we describe the vocabulary associated with the invention, before presenting the drawbacks of the prior art, and then finally showing in greater detail how the invention remedies them.

[0121] As used herein, “metadynamics based method” can refer to any collective-variable-based enhanced-sampling technique of the MetaDynamics family (including well-tempered MetaDynamics and its OPES / OPES-Explore variants) in which a history-dependent biasing energy term is built on-the-fly so as to drive sampling toward a prescribed target distribution (e.g., a well-tempered distribution) and facilitate transitions across free-energy barriers. As used herein, “adaptive biasing method” can refer to any enhanced-sampling method that updates, during simulation, a biasing force and / or biasing energy term as a function of accumulated sampling statistics along one or more collective variables.

[0122] As used herein, a “biasing force” is a force that is added to an original physical Hamiltonian. In some embodiments, a biasing force is added on top of an alchemical Hamiltonian for example to guide sampling along A.

[0123] As used herein, a “bias potential” is an artificial potential energy term added to the molecular system to alter a sampling distribution of a chosen variable, such as the alchemical parameter A or another CV. In some embodiments, the bias potential is updated on the fly during simulation.

[0124] As used herein, “collective variable (CV)” can refer to a reduced-dimensional coordinate (scalar or vector) defined as a function of the microscopic degrees of freedom of a molecular system, selected so that biasing forces, biasing energy terms, or probability reweighting can be applied along it to facilitate sampling of the configuration space. In particular, the alchemical coupling parameter A can be treated as a collective variable.

[0125] As used herein, “Alchemical parameter (A)” or “Alchemical coupling parameter (A)” can refer to a variable which can control interpolation between two end-state Hamiltonians (e.g., a ligand fully non-interacting versus fully interacting). In some embodiments, A is between 0 and 1. In some embodiments, reflective boundary conditions are defined for 'A=0' and 'A=1' to ensure A remains within [0,1] during simulation.As used herein, a “free-energy surface” (FES) can refer to a free-energy function expressed as a function of one or more collective variables of the molecular system. In some embodiments, the FES is computed from simulation data collected while sampling the CV(s) under an adaptive bias, and is obtained by thermodynamic integration of energy derivatives or mean forces along the CV(s). In particular, when the CV comprises the alchemical coupling parameter A, the free-energy surface can correspond to the alchemical free energy as a function of A.

[0126] As used herein, “flattened distribution” can refer to a sampling distribution along a CV (e.g., A) in which free-energy barriers along that coordinate are reduced, such that the sampled distribution approaches a substantially uniform distribution along that coordinate. In some embodiments, ABF applies biasing forces along the transition coordinate to flatten the sampled free-energy landscape toward a uniform distribution, or substantially uniform distribution.

[0127] As used herein, “sampling” can refer to how a molecular simulation explores configuration space and, when present, an alchemical parameter space. In some embodiments, sampling is performed in a single simulation over a continuous range of a collective variable (e.g., A), rather than by splitting that range into discrete windows, by treating A as an extra degree of freedom (expanded-ensemble / lambda-dynamics scheme).

[0128] As used herein, “molecular system” refers to a set of atoms that constitutes the entire simulation object, including any solute, solvent, gas, or other components that may be present. For example a molecular system can be selected among : ligands, small molecular fragments (inorganic or organic), peptides or short proteins, nucleic-acid fragments, cofactors, metals, or metal-coordinated complexes. In some embodiments, the molecular system comprises a ligand and a receptor and is simulated under periodic boundary conditions.

[0129] As used herein, a “target probability distribution” can refer to a prescribed marginal probability distribution along one or more collective variables that the biased simulation is configured to approach (i.e., the bias potential and / or biasing force is updated so that the sampled CV distribution approaches the target distribution). In some embodiments, the target probability distribution is a well-tempered distribution.

[0130] As used herein, “well-tempered distribution” can refer to a target probability distribution used in probability-based enhanced sampling methods that is associated with reduced free-energy barriers while preserving thermodynamic relevance. In some embodiments, OPES-Explore broadens sampling toward a target probability distribution known as the well-tempered distribution. In some embodiments, the well-tempered distribution is characterized by a “bias factor” parameter that controls the extent of smoothing / broadening of the sampled distribution.As used herein, an “equilibrium (unbiased) probability distribution” of a collective variable can refer to the marginal probability distribution that would be obtained for that collective variable in the absence of the applied bias (i.e., under the underlying physical or alchemical Hamiltonian without the added bias potential). In some embodiments, the adaptive biasing force is computed from an on-the-fly estimate of such an equilibrium (unbiased) distribution.

[0131] As used herein, a “current biased probability distribution” can refer to the marginal probability distribution of a collective variable observed under the bias applied during the simulation (e.g., under the current bias potential / adaptive biasing force). In some embodiments, the adaptive biasing force is computed from a biasing potential issued from an on-the-fly estimate of such current biased probability distribution, rather than directly from an equilibrium target distribution.

[0132] As used herein, an “on-the-fly estimate” can refer to an estimate (e.g., of a probability distribution, mean force, or bias potential) that is updated during execution of the simulation as additional sampling data is accumulated, rather than being determined only after completion of the simulation (post-processing).

[0133] As used herein, “thermodynamic integration” (Tl) can refer to computing a free-energy surface or free-energy difference by integrating, over a range of a collective variable (preferably A), energy derivatives of the molecular system’s potential energy with respect to that collective variable. In some embodiments, Tl comprises continuously integrating such ensemble-averaged derivatives over A to obtain an alchemical free energy.

[0134] As used herein, an “energy derivative” (or “potential energy derivative”) with respect to a collective variable can refer to a derivative of the molecular system’s potential energy with respect to that collective variable, computed during the simulation from collected information. In some embodiments, the energy derivative is 3U / 3A and is used within thermodynamic integration.

[0135] As used herein, a “mean force” along a collective variable can refer to an ensemble-averaged force associated with that collective variable (e.g., an averaged derivative of free energy with respect to the collective variable, such as along A). In some embodiments, mean forces exerted along A are numerically integrated to obtain a free-energy difference between a bound state and an unbound state.

[0136] As used herein, a “free-energy gradient” can refer to a derivative of a free-energy surface with respect to a collective variable (e.g., dF / dA). In some embodiments, multiple-walker simulation data are aggregated to compute such a free-energy gradient for each value (or region) of A, and the aggregated mean force / gradient is integrated to obtain the alchemical free energy.As used herein, “Gaussian kernel” can refer to Gaussian functions used to construct, represent, and / or update an on-the-fly bias potential. It can be used to build and / or represent (i) a bias potential and / or (ii) an estimated probability distribution of the collective variable(s). In some embodiments, OPES employs Gaussian kernels to adaptively construct a bias potential that guides exploration while preserving thermodynamic relevance.

[0137] As used herein, a “Gaussian kernel bandwidth” (or “Gaussian kernel width”) can refer to a parameter that determines the spread / width of Gaussian kernels used in the kernel-based representation (e.g., for estimating a distribution and / or constructing a bias potential). In some embodiments, the Gaussian kernel width is determined using an adaptive sigma. As used herein, “adaptive sigma” can refer to an on-the-fly updated kernel-width selection mechanism used to determine Gaussian kernel width during the simulation, optionally varying as sampling proceeds.

[0138] As used herein, a “kernel deposition frequency” (or “kernel deposition pace”) can refer to how often Gaussian kernels are deposited / added / updated during a simulation, expressed as a number of simulation steps and / or a deposition time interval. In some embodiments, kernels are deposited at a fixed step interval (e.g., 300 steps).

[0139] As used herein, a “barrier parameter” (AE) can refer to a parameter used in OPES / OPES-Explore that limits and / or regularizes the applied bias so as to (i) allow efficient transitions between basins and (ii) prevent access to irrelevant high-energy states. In some embodiments, the barrier parameter is selected as a bias threshold sufficient to push the system out of local minima, and need not equal the free-energy difference associated with the simulation.

[0140] As used herein, a “weighted kernel density estimation” can refer to an estimation procedure for an underlying probability distribution of a collective variable based on simulation samples, wherein each sample is associated with a weight related to the current bias so as to approximate an unbiased probability distribution.

[0141] As used herein, “kernel compression” (or “compressing the kernel representation”) can refer to limiting computational growth of a kernel-based distribution representation by merging kernels corresponding to previously sampled regions, thereby refining a distribution estimate over time while maintaining bounded memory usage.

[0142] As used herein, an “adaptive kernel bandwidth” or “adaptive kernel width’can refer to a kernel width selection scheme in which the kernel bandwidth decreases as the number of effective samples increases, providing a coarser estimate early in a simulation and progressively higher resolution as additional data accumulate, thereby improving bias accuracy and stability without requiring manual parameter tuning.

[0143] As used herein, a “number of effective samples” can refer to an estimate of the amount of statistically useful sampling data contributing to the distribution estimate (including in aweighted setting), used to control resolution (e.g., kernel bandwidth) of the probability estimate and / or bias construction as sampling proceeds.

[0144] As used herein, “multiple-walker sampling instances” can refer to a plurality of concurrently executed simulation replicas (“walkers”) used within a multiple-walker framework, each walker generating its own trajectory while collectively contributing to the sampling. Preferably they are configured to share running estimates used by a collectivevariable biasing method. “Multiple-walker sampling” thus refers to multiple concurrent sampling instances that may exchange or share mean-force and / or bias information.

[0145] As used herein, a “fully coupled state” can refer to an alchemical end state in which the ligand is fully interacting with its environment (e.g., protein and / or solvent), and a “decoupled state” can refer to an alchemical end state in which such interactions are scaled down (optionally to zero), with A interpolating between the two.

[0146] As used herein, “interpolation of polarization between end states” can refer to computing potential derivatives with respect to A for polarizable force fields by interpolating polarization contributions between alchemical end states such that the associated derivative can be obtained as a difference between end-state polarization energies.

[0147] As used herein, “Distance-to-Bound-Configuration (DBC)” is a restraint variable used to define a bound state in binding free energy calculations. In some embodiments, DBC is defined as an root mean square deviation of selected ligand atoms with alignment to selected receptor binding-site atoms.

[0148] As used herein, “DBC restraint” can refer to a flat-bottom harmonic restraint applied on the Distance-to-Bound-Configuration of each molecular system using two independent, A-dependent target distances to control the restraint strength throughout the alchemical transformation.

[0149] As used herein, a “DBC cutoff” can refer to a threshold value for the Distance-to-Bound-Configuration variable selected from a monitored DBC distribution (e.g., monitored during plain molecular dynamics), used as a cutoff for applying a DBC restraint. In some embodiments, the DBC cutoff is selected as a DBC value within a 95% interval of the monitored DBC distribution.

[0150] As used herein, “plain molecular dynamics” (plain MD) can refer to an unbiased molecular dynamics simulation of a system performed without applying the enhanced-sampling biasing method (e.g., without applying the Lambda-ABF-OPES bias along A), and used to monitor distributions (e.g., DBC) and / or atom fluctuations (e.g., RMSF) for selecting restraint parameters.

[0151] As used herein, a “flat-bottomed harmonic restraint” can refer to a restraint potential that is substantially inactive (flat) below a defined cutoff of a restraint variable and becomesharmonic above that cutoff. In some embodiments, a flat-bottomed harmonic restraint is applied to the DBC above the DBC cutoff.

[0152] As used herein, “root mean square fluctuation (RMSF)” can refer to a statistical measure of positional fluctuations of atoms over a molecular dynamics trajectory, used to select atoms for defining a restraint variable (e.g., selecting ligand heavy atoms and / or receptor anchor atoms for DBC based on an RMSF threshold).

[0153] As described hereafter, the inventors developed a new method for computing a free-energy surface and, in particular, an alchemical free-energy difference (e.g., an alchemical binding free energy) of a molecular system in a molecular dynamics simulation.

[0154] This method involves combining an Adaptive Biasing Force (ABF) scheme with an OPES-type bias to enhance sampling along one or several collective variables (CVs). Preferably, this method involves treating an alchemical parameter A as a collective variable (CV) and combining an Adaptive Biasing Force (ABF) scheme with an OPES-type bias to enhance sampling along the alchemical parameter A. As illustrated in examples, OPES-Explore alone struggles to compute the binding free energies of interest and final AG values varied between runs. However, the particular hybrid method of the invention lead to reduced compute time and improved convergence and robust simulation workflow without postprocessing.

[0155] In particular, the inventors developed the use of a bias potential expressed as a function of a target probability distribution (preferably a well-tempered distribution) and constructed / updated using Gaussian kernels, in combination with an adaptive biasing force acting on the collective variables (CVs), preferably the alchemical parameter A. Such features allow for improved exploration of the alchemical coordinate, mitigation of kinetic trapping, and faster convergence of free-energy estimates while maintaining accuracy.

[0156] As it will be described, the alchemical free-energy difference can be computed by thermodynamic integration (Tl) of energy derivatives with respect to collective variables (CVs), preferably A, collected during the biased sampling.

[0157] Also the invention preferably relies on a portable computer implementation (e.g., via the CoIvars library) and / or multiple-walker sampling instances that exchange or share biasing information.

[0158] Advantageously, this approach involves continuously sampling the alchemical parameter A as a dynamical variable over its domain (e.g., 0 < A < 1 , optionally with reflective boundary conditions) rather than using multiple fixed-A windows. This facilitates on-the-fly monitoring of convergence (e.g., by evaluating partial integrals of the mean force with respect to the alchemical parameter A over cumulative simulation time) and can reduce the need for postprocessing.This approach is particularly advantageous for computing binding affinity / absolute binding free energy of a ligand to a target, and more generally for computing free-energy differences and profiles in molecular systems (including, for example, solvation free energies, conformational free-energy differences, etc.).

[0159] Hence, according to a first aspect, the invention relates to a computer-implemented method 100 for computing a free-energy surface (and / or an alchemical binding free-energy surface / alchemical free-energy difference) in a molecular system. In particular, the invention relates to a computer-implemented method for computing an alchemical free-energy difference between two end-state Hamiltonians of a molecular system by molecular dynamics.

[0160] This method can be used for computing a binding free energy of a ligand to a target / an absolute binding free energy / a binding affinity, and / or other free-energy computations performed by molecular dynamics simulation. The computer-implemented method preferably comprises the use of an operating computer system comprising one or more processors and one or more non-transitory computer-readable memories storing instructions for carrying out the method.

[0161] The method can comprise the following steps:

[0162] - Obtaining (and optionally defining) one or several collective variable(s) describing the molecular system, wherein said CV(s) comprise an alchemical parameter A;

[0163] - Applying, during a molecular dynamics simulation, to the alchemical parameter A:

[0164] o an adaptive biasing force (ABF) updated on the fly so as to drive sampling of A toward a flattened (e.g., substantially uniform) distribution, and o a bias potential as an additional potential energy term updated on the fly and defined as a function of a target probability distribution (preferably a well- tempered distribution) of the alchemical parameter A, the bias potential being constructed using Gaussian kernels.

[0165] - Computing the alchemical free-energy difference (or free-energy surface) by thermodynamic integration of energy derivatives with respect to A obtained / collected during step 2.

[0166] Optionally, the method further comprises running multiple walkers and exchanging / sharing biasing information between walkers.

[0167] The figure 1 illustrates a particular embodiment of the method according to the invention. As shown in figure 1, a method according to the invention can comprise a step of obtaining 110 at least one collective variable describing the molecular system.

[0168] This step 110 can be designed to define a reduced set of coordinates that capture the slow degrees of freedom governing the thermodynamic transformation of interest, and in particularto introduce an alchemical coupling parameter enabling a continuous transformation between two predefined thermodynamic states of the molecular system. For example, the at least one collective variable comprising an alchemical parameter A that interpolates between the two end-state Hamiltonians (e.g., a ligand fully non-interacting versus fully interacting).

[0169] In some embodiments, there is only one collective variable. In some embodiments, there are several collective variables. In some embodiments, A is between 0 and 1. In some embodiments, reflective boundary conditions are defined for A=0 and A=1 to ensure that A remains within [0,1] during the simulation.

[0170] In particular, the at least one collective variable comprises an alchemical parameter A that interpolates between the two end state Hamiltonians, such that a first Hamiltonian corresponding to A = 0 defines a first physical state of the system and a second Hamiltonian corresponding to A = 1 defines a second physical state of the system, the intermediate values of A defining a continuous alchemical pathway between said states.

[0171] Preferably, the alchemical parameter A is defined as a scalar variable bounded within a predefined interval, in particular [0,1], and is associated with a A-dependent Hamiltonian constructed by scaling selected interaction terms of the molecular system, including nonbonded electrostatic and / or van der Waals interactions of a ligand or molecular fragment relative to its environment.

[0172] In some embodiments, the step of obtaining 110 the collective variable further comprises:

[0173] identifying a subset of atoms forming an alchemically transformed region; defining the dependence of the system Hamiltonian on A through scaling functions applied to interaction parameters of said subset; and

[0174] initializing A as a dynamical variable suitable for continuous propagation during a molecular dynamics simulation.

[0175] Advantageously, defining A as a continuous collective variable enables sampling of the alchemical space without discretization into predefined A windows, thereby facilitating adaptive biasing and thermodynamic integration within a single simulation framework.

[0176] As shown in figure 1, a method according to the invention can comprise a step of updating 120 on the fly and applying to the collective variables (CVs), preferably the alchemical parameter A a biasing scheme.

[0177] This step can be designed to progressively modify the effective free-energy landscape along the alchemical coordinate so as to enhance transitions between intermediate alchemical states, reduce kinetic trapping in A-space, and ensure efficient and uniform exploration of the alchemical pathway during the molecular dynamics simulation.

[0178] A biasing scheme according to the invention comprises the use of an adaptive biasing force 121 and a bias potential 122.The adaptive biasing force is preferably configured to drive sampling of A toward a flattened distribution. This is preferably done by estimating on the fly the mean force acting along A and applying a compensating force so as to reduce free-energy gradients associated with the alchemical transformation.

[0179] In some embodiments, the adaptive biasing force is computed from an on-the-fly estimate of an equilibrium (unbiased) target probability distribution. In some embodiments, the adaptive biasing force is computed from a biasing potential issued from an on-the-fly estimate of a current biased probability distribution (wherein a probability density observed under the bias is used to update the bias rather than a target (equilibrium) distribution). In some embodiments, the adaptive biasing force is computed from a biasing potential issued from an on-the-fly estimate of a well-tempered distribution being sampled

[0180] The bias potential is preferably added as an additional potential energy term to the Hamiltonian of the molecular system, such that the total potential energy becomes the sum of the physical interaction energy and the bias contribution acting on A. Preferably, the bias potential is defined as a function of a target probability distribution of A. Alternatively it can be defined as a function of an on-the-fly estimate of the sampled probability distribution of A, or as a function of a reweighted estimate of the underlying free-energy profile along A. In some embodiments, step b) employs weighted kernel density estimation to estimate an underlying probability distribution of the biased variable (e.g., A) and to adjust the bias potential so that a simulated distribution of the biased variable and / or CV(s) approaches a target probability distribution

[0181] Advantageously, the target probability distribution is a well-tempered distribution, configured to reduce effective free-energy barriers while preserving thermodynamic consistency and enabling recovery of unbiased statistics. Alternatively, the target probability distribution can be a uniform distribution over a predefined A interval, or another predefined non-uniform distribution selected to control exploration of specific alchemical regions.

[0182] As describe in examples, the bias potential is advantageously constructed using Gaussian kernels. Preferably, the parameters of said kernels being updated on the fly so as to progressively approximate the target probability distribution and to regulate the explorationconvergence balance of the enhanced sampling scheme. In some embodiments, Gaussian kernels are defined with parameters comprising at least one from: a Gaussian kernel bandwidth (sigma), a bias factor, and a barrier parameter 'AE' (preferably the barrier parameter). In some embodiments, an adaptive sigma is used to determine Gaussian kernel width. In some embodiments, a frequency for kernel deposition is set as at least one of: a number of simulation steps and a deposition time interval (see §5.5). In one non-limiting example configuration, a deposition frequency is set to 300 steps, corresponding to about 9ps when a time step is 3 fs. In one non-limiting example configuration, a bias threshold (barrier parameter) is set to at least 5 kcal / mol; in another non-limiting example configuration (e.g., a weak-binder case), a bias threshold of 2 kcal / mol is used.

[0183] As a detailed description of the implementation of a method of the invention, in some embodiments, step b) comprises initializing initial conditions for the molecular system, including atomic coordinates, velocities, and an initial value of the biased variable (e.g., A), as well as initializing any bias-state variables (e.g., ABF accumulators, kernel lists, and / or distribution estimates) used to update the bias on-the-fly. In some embodiments, preparing the simulation includes generating a solvated and neutralized system prior to biasing on A. In one non-limiting example configuration, a molecular system comprises on the order of tens of thousands of atoms (e.g., about 31,000-42,000 atoms), optionally includes a physiological salt concentration (e.g., 0.15 M NaCI), and is placed in a periodic water box (e.g., a cubic box of about 70-80 A) with a padding distance (e.g., at least about 12 A between a solute and a box edge).

[0184] In some embodiments, molecular system coordinates (and, in some embodiments, velocities) are propagated by molecular dynamics and / or lambda-dynamics to generate a trajectory and sample configuration space and biased-variable space. In some embodiments, A is treated as a dynamical variable whose value changes during the simulation and is biased as described herein.

[0185] In some embodiments, information is collected about the molecular system’s energy to support ABF updates and Tl, including energy derivatives of the molecular system’s potential energy with respect to A and / or mean forces exerted along A. In some embodiments, such quantities are accumulated and / or stored as a function of A (e.g., in bins or by kernel-based representations) for subsequent integration and / or monitoring.

[0186] In some embodiments, an adaptive biasing force (ABF) scheme adaptively computes a derivative of a free energy associated with the parameter A using a thermodynamic integration (Tl) formula, and the estimate is applied as a force directly to the simulated system to guide the dynamics. In some embodiments, as a result of applying ABF, sampling of A converges toward a uniform distribution, facilitating the crossing of free-energy barriers and driving sampling of A toward a flattened distribution.

[0187] In some embodiments, the bias is not predetermined; instead, the simulation builds an adaptive bias on-the-fly by measuring the system’s response. In some embodiments, the adaptive bias counter-balances a free-energy profile along A, making A values equally probable, and unbiased free-energy surfaces can be recovered after the simulation by removing the known bias.

[0188] In some embodiments, the biasing scheme includes use of a bias potential that is an artificial potential energy term added to the system to alter the sampling distribution of a chosenvariable (e.g., A or another CV). In some embodiments, the bias potential is expressed as a function of a target probability distribution of the biased variable (e.g., A). In some embodiments, the target probability distribution is a well-tempered distribution. In some embodiments, the bias potential is constructed using Gaussian kernels and is updated on-the-fly.

[0189] In some embodiments, a kernel-based enhanced sampling approach (e.g., OPES-Explore as a non-limiting example) broadens the sampling of a system toward a target probability distribution known as a well-tempered distribution by adaptively constructing a bias potential using Gaussian kernels. In some embodiments, parameters for such a kernel-based approach include a Gaussian kernel width (sigma) and a bias factor, and such parameters are used to control the sampling scope and stability. In some embodiments, a barrier parameter AE is employed to allow transitions between basins while preventing access to irrelevant high-energy states.

[0190] In some embodiments, a combined biasing comprises the ABF and the bias potential, and the combined biasing is applied during the simulation to drive sampling of the biased variable (e.g., A) toward a flattened distribution. In some embodiments, ABF reduces free-energy barriers, while the kernel-based bias potential incorporates a history-dependent potential term that pushes the system out of local minima.

[0191] In some embodiments, information collected during step b) is stored and / or aggregated as a function of the biased variable (e.g., A) for Tl, including mean-force and / or energy-derivative information. In some embodiments, Tl data is accumulated per walker and / or aggregated across walkers for subsequent integration.

[0192] In some embodiments, step b) is applied concurrently by multiple-walker sampling instances. In particular, step b) can be applied concurrently by multiple-walker sampling instances that exchange, share, transmit, and / or aggregate mean-force and / or bias information. In some embodiments, multiple-walker sampling instances exchange or share adaptive biasing force data and / or kernel-based bias information so as to collectively explore different regions of phase space. In some embodiments, each multiple-walker sampling instance periodically transmits locally accumulated biasing data to a central data structure or to all other multiplewalker sampling instances. In one non-limiting example configuration, at least four walkers are employed, preferably four walkers are employed.

[0193] In some embodiments, free-energy is monitored in real time by evaluating partial integrals of the mean force with respect to A and analyzing resulting estimates over cumulative simulation time.

[0194] In some embodiments, sampling in the context of the invention refers to how the molecular simulation explores both configuration space and biased-variable space. In some embodiments, the method is designed to continuously sample the full range of the biasedvariable (e.g., 0 < A < 1) in a single simulation, rather than breaking the range into discrete windows.

[0195] As shown in figure 1, a method according to the invention can comprise a step of computing 130 the alchemical free-energy difference.

[0196] This step can be designed to determine, from the enhanced sampling performed along the collective variables (CVs), preferably the alchemical parameter A, the free energy variation between two states. The two states are for example the two end-state Hamiltonians corresponding to A = 0 and A = 1 , thereby quantifying the thermodynamic difference between the associated physical states of the molecular system.

[0197] In particular this step can imply computing the thermodynamic integration of energy derivatives calculated from the collected information with respect to A obtained during the step b). More specifically, the method comprises estimating, during the biased molecular dynamics simulation, the ensemble-averaged derivative of the system Hamiltonian with respect to A, and integrating said derivative over the predefined A interval.

[0198] Advantageously, the thermodynamic integration is performed continuously within a single simulation in which A is propagated as a dynamical variable, such that the free energy difference is obtained without requiring discretization into multiple fixed-A windows and without post-processing through separate free energy estimators.

[0199] For example, the alchemical free energy difference AG can be computed as the integral, over A from 0 to 1 , of the mean force (5H / <5A)A accumulated during the biased sampling, the integral being evaluated numerically using the on-the-fly averaged derivatives collected during the simulation.

[0200] In some embodiments, in step c), thermodynamic integration (Tl) is performed to compute a free-energy surface and / or an alchemical free-energy difference between end states, and an output is generated. In some embodiments, Tl comprises integrating, with respect to A, energy derivatives of the molecular system’s potential energy with respect to A and / or mean forces exerted along A, over a range of A values sampled during step b).

[0201] In some embodiments, Tl comprises continuously integrating an ensemble-averaged derivative of potential energy with respect to A over a range of A values sampled continuously in a single A-dynamics simulation, without breaking the A range into discrete fixed-A windows. In some embodiments, step (c) is carried out during step (b) so as to provide an on-the-fly estimate of an alchemical free-energy difference during the simulation.

[0202] In some embodiments, step (c) comprises computing a free-energy profile as a function of A by integrating mean forces with respect to A, and obtaining an alchemical free-energy difference between end-state Hamiltonians from the integrated profile.In some embodiments, tree-energy is monitored in real time by evaluating partial integrals of the mean force with respect to A and analyzing resulting estimates over cumulative simulation time.

[0203] Preferably, when multiple-walker sampling is used, mean-force and / or energy-derivative information collected during step b) is aggregated across walkers and integrated to obtain an alchemical tree-energy estimate.

[0204] An output can comprise a tree-energy surface and / or a tree-energy difference, including, in some embodiments, an alchemical binding free energy. In some embodiments, the output is generated during the simulation (e.g., on-the-fly) and / or after completion of the simulation. In some embodiments, a Distance-to-Bound-Configuration (DBC) coordinate is used to restrain a ligand during alchemical simulations, including, in some embodiments, during equilibration. In some embodiments, DBC is defined as a root mean square deviation (RMSD) of selected ligand atoms for each frame, with alignment performed relative to selected atoms of a receptor binding site. This can capture positional, orientational, and conformational deviations of the ligand in a single collective variable.

[0205] In some embodiments, selecting ligand atoms used for the DBC variable comprises monitoring a root mean square fluctuation (RMSF) of ligand heavy atoms over a plain molecular dynamics simulation (e.g., a simulation of at least 100 ns) and selecting atoms based on RMSF. In one non-limiting example configuration, ligand atoms with an RMSF of less than 0.6 A are selected for tight binders, while ligand atoms with an RMSF between 0.7 A and 0.8 A are selected for weak binders.

[0206] In some embodiments, selecting receptor binding-site atoms comprises selecting protein atoms within a distance of the ligand (e.g., selecting Co atoms within 6 A of the ligand). In one non-limiting example configuration, selected receptor-side Co atoms have RMSF values of about 0.6 A. In some embodiments, during plain molecular dynamics of a ligand, DBC is monitored, and a DBC value within a 95% interval of a DBC distribution is selected as a DBC cutoff for a restraint to reduce restraint-induced artifacts.

[0207] In some embodiments, a flat-bottomed harmonic restraint is applied to the DBC above the DBC cutoff. In one non-limiting example configuration, a flat-bottomed harmonic restraint is applied above the DBC cutoff with a force constant of 100 kcal / mol / A2.

[0208] In some embodiments, a system is provided for performing the described method(s).

[0209] Hence, in another aspect, the invention relates to a system configured to implement the method according to the invention. In some embodiments, the system comprises one or more processors and one or more non-transitory computer-readable memories storing instructions which, when executed, cause the processor(s) to perform any of the methods described herein.In some embodiments, the system comprises a molecular dynamics engine and a CV / biasing module configured to implement the on-the-fly ABF and kernel-based bias potential described herein.

[0210] In one non-limiting example implementation, molecular dynamics simulations are performed using Tinker-HP (e.g., Tinker-HP Version 1.2) including the CoIvars library for Lambda-ABF-OPES simulations, and a kernel-based enhanced sampling technique (e.g., OPES-Explore) is provided by the CoIvars library (e.g., a CoIvars version dated 2024-11-18). In some embodiments, such an integrated environment provides convergence estimation without requiring post-processing and supports compatibility with other CV-based methods.

[0211] In one non-limiting example implementation, plain OPES simulations are performed using Tinker-HP and Plumed2.

[0212] In some embodiments, the system is configured for multiple-walker operation in which multiple computing instances concurrently sample A and exchange, transmit, and / or aggregate mean-force and / or bias information.

[0213] In some embodiments, a non-transitory computer-readable medium stores instructions which, when executed by a processor, cause the processor to perform any of the described methods.

[0214] The described methods and systems are applicable to molecular simulation workflows including, for example, drug discovery and design, protein-ligand binding free energy estimation, and other applications where free-energy surfaces and differences are computed. Non-limiting examples include computing: a binding free energy of a ligand to a target; an absolute binding free energy of a ligand; a binding affinity between a ligand and a target; a binding affinity between a ligand and a receptor; a free-energy profile for protein folding; free-energy differences between molecular conformational states; activation free energies; reaction free energies; solvation free energies; partition coefficients; dissociation constants; association constants; relative stabilities of molecular conformers; free-energy differences in enzyme-substrate complexes; free energies for conformational transitions; free energies for phase transitions; free energies of adsorption; free energies for ligand-induced conformational changes; free energies in molecular docking simulations; free energies for applications in drug design; free energies for stability prediction of molecular systems; free energies in polymer systems; free energies in solid-state systems; free energies in solvated systems; free energies for protein-ligand interactions; free energies for nucleic acid-ligand interactions; free energies in molecular dynamics simulations; free energies for conformational sampling; and / or free energies for estimating reaction rates.EXAMPLE

[0215] The invention is further described in detail by reference to the following experimental examples. These examples are provided for purposes of illustration only and are not intended to be limiting unless otherwise specified. Thus, the invention should in no way be construed as being limited to the following examples, but rather, should be construed to encompass any and all variations which become evident as a result of the teaching provided herein.

[0216] Without further description, it is believed that one of ordinary skill in the art can, using the preceding description and the following illustrative examples, make and utilize the quantum circuit and ansatz of the present invention and practice the claimed methods. The following working examples therefore, specifically point out the preferred embodiments of the present invention, and are not to be construed as limiting in any way the remainder of the disclosure. Numerical values are provided as example settings observed / used in the described simulations and are not limiting unless explicitly recited in the claims.

[0217] Methods

[0218] MD Simulations

[0219] We performed all simulations using the molecular dynamics code Tinker-HP Version 1.2, which includes the CoIvars library (Fiorin et al., 2013 ) for Lambda-ABF-OPES simulations. For the protein, water, and ion parameters, the polarizable AMOEBA

[0220] force field was applied. Ligand parameterization was performed using the Poltype package. The systems, containing between -31,000 and -42,000 atoms on average (varying with ligand size), were neutralized with 0.15 M physiological NaCI concentration. A cubic water box with dimensions of 70-80 A was used, ensuring a minimum distance of 12 A between the ligand and the box edge, with its size adjusted according to the ligand.

[0221] The systems were first minimized using the Tinker-HP minimize program, followed by a two-step heating process from 200 K to 298 K under the NVT ensemble. In the first step, restraints were applied to all protein-ligand and X-ray water molecule atoms. In the second step, restraints were applied to heavy backbone atoms, using a force constant of 10 kcal / mol / A2in both steps. Each stage included 4 ns of MD simulations, conducted using the RESPA integrator with a time step of 2.0 fs.

[0222] In the third step, the three-level multiple timestep BAOAB-RESPA1 integrator with an outer timestep of 10 fs (and an intermediate of 3.33 fs and a shorter timestep of 1 fs) was employed under the NVT ensemble for 5 ns, with restraints applied only to Ccratoms using a force constant of 1 kcal / mol / A2. This was followed by equilibration in the NPT ensemble at 1 atm for 5 ns, maintaining the same restraints as the previous step. Finally, for each ligand, a production run of at least 100 ns was performed under the NPT ensemble using the BAOAB-RESPA1 Langevin integrator with a 10 fs time step. During the production phase, a restraint with a force constant of 1 kcal / mol / A2was applied to the Ccr atoms of the protein’s floppy tail to prevent large conformational changes.

[0223] Non-Langevin temperature control was ensured by using the Bussi thermostat and pressure with the Berendsen barostat Van der Waals interactions employed a 12 A cutoff, while electrostatic interactions were treated using the Particle Mesh Ewald (PME) method with a real-space cutoff of 7 A. Induced dipoles were calculated with a Preconditioned Conjugate Gradient (PCG) solver, with a convergence tolerance of 1 x 10“5Debye.

[0224] Analytical lambda derivatives for the polarization energy using Particle Mesh Ewald A simple interpolation of polarization between the end states reads:

[0225]

[0226] so that:

[0227] dE

[0228] (r, A) = Epolr, 1) - Epolr, 0)

[0229] uA

[0230] which requires two resolutions of the polarization equations per timestep as stated in the main text. Alternatively, let’s consider the polarization energy associated to a scaling of the polarizabilities and the permanent multipoles of the "alchemical" part of the system, the complete system being made of N atoms:

[0231]

[0232] where T is the (3N,3N) polarization matrix and E the 3N vector of the permanent electric fields on the polarizable sites. The Hellman-Feynman theorem yields:

[0233]

[0234] The derivative contribution due to the T is trivial to compute because it is only associated to its diagonal part «“1(A), with « collecting the diagonal polarizability tensors. In the context of periodic boundary conditions computed with Particle MeshEwald, the second one can be separated in 3 given the various component of — = uA

[0235]

[0236] The first two terms can be directly computed. If only fixed charge are used (and no permanent multipoles) then the last one can be reformulated as:

[0237]

[0238] where Q A) is the vector (of size N) containing the permanent charges and VreciPthe vector of same size containing the reciprocal potential due to the induced dipoles. This potential is readily available to compute the permanent polarization energy and dQW is trivial to compute. The same reasoning is naturally extended to permanent dA

[0239] multipoles of higher order by involving derivatives of VredP(ji) which are always available to compute the permanent polarization energy.

[0240] DBC Restraint

[0241] In alchemical methods, absolute binding free energy calculations critically depend on a precise definition of the bound state and well-designed ligand restraints (both translational and, optionally, orientational) to ensure rapid convergence. In this implementation, we employed the DBC coordinate to restrain the ligand during Lambda-ABF-OPES simulations. It is defined as the RMSD of some ligand atoms (LA) for each frame, with the alignment performed relative to some atoms of the receptor’s binding (RB) site. This approach captures positional, orientational, and also conformational deviations of the ligand in a single collective variable. To select the LA, we monitored the root mean square fluctuation (RMSF) of the ligand’s heavy atoms over at least 100 ns of standard MD simulations. Atoms with an RMSF of less than 0.6 A were selected for tight binders (Ligands 1 -9), while those with an RMSF between 0.7 and 0.8 A were chosen for weak binders (Ligands 10-11 ). For the RB selection, we used the Ccr atoms of the protein within 6 A of the ligand, with an RMSF of approximately 0.6 A. See Fig. 8 for the definition of DBC for the protein and each ligand. During plain MD of each ligand, the DBC was monitored using the CoIvars library, and the DBC value within the 95% interval of the distribution was selected as the DBC cutoff for the restraint in Lambda-ABF-OPES. This selection ensures thatduring the alchemical simulation, there is no biased artifact from the DBC restraint.

[0242] The DBC cutoff for each ligand is listed in Table 1. A flat-bottomed harmonic restraint was then applied to the DBC with a force constant of 100 kcal / mol / A2above the DBC cutoff.

[0243] Table 1

[0244] Compound DBCCTii{7(A) DBC (kcal / mol) Harmonic (kcal / mol) L1-BRD4 0.91 4.52 ± 0.07 2.33

[0245] L2-BRD4 0.79 5.00 ± 0.05 2.51

[0246] L3-BRD4 0.80 5.14 ± 0.03 2.43

[0247] L4-BRD4 0.91 3.77 ± 0.05 2.33

[0248] L5-BRD4 0.87 4.64 ± 0.00 2.37

[0249] L6-BRD4 0.86 5.47 ± 0.02 2.38

[0250] L7-BRD4 0.84 5.20 ± 0.15 2.40

[0251] L8-BRD4 0.88 4.88 ± 0.02 2.36

[0252] L9-BRD4 0.84 5.37 ± 0.00 2.39 L10-BRD4 0.96 4.19 ± 0.04 2.28 L11-BRD4 0.97 5.15 ± 0.04 2.28 Benzamidine-Trypsin 0.81 4.68 ± 0.03 2.42

[0253] Restraining specific regions of a protein, such as loop regions or flexible side chains, is preferred to prevent them from adopting unrealistic conformations that could distort the calculated free energy. This strategy is especially critical when a binding event stabilizes a particular protein conformation that needs to be accurately maintained throughout simulations.

[0254] We applied a positional restraint with a mild force constant of 1 kcal / mol / A2on the Co atoms of the flexible tail of BRD4 during plain MD simulations. For Lambda-ABF-OPES, positional restraints were applied to all Co atoms as well as the heavy backbone atoms within 6 A of the ligand and the heavy atoms of Asn140 with a force constant of 2 kcal / mol / A2. The harmonic restraint was activated only when the displacement exceeded 0.5 A, ensuring that the normal fluctuations of the atoms were preserved while preventing excessive deviations during the simulation.

[0255] Atoms Involved in the Definition of DBC for Restraint and DBC Correction Atoms involved in the DBC restraint for the Lambda-ABF-OPES simulations in the complex phase are shown in Fig. 8. The Ca atoms within 6 A of the ligand, with RMSF values less than 0.6 A, used as anchor atoms for the DBC, are highlighted in yellow, and the corresponding residue numbers are shown. Notably, for L2, two additional atoms are included in the definition of the DBC anchor atoms, with indexnumbers 6 and 7 (Fig. 8(a)). For each ligand, the atoms selected for the DBC are highlighted in orange (Fig. 8(b)). These correspond to atoms with RMDF values less than 0.6 A for L1 to L9 and less than 0.8 A for L10 and L11 , which are the tightest binders in the list and relax within the pocket over 100 ns of plain MD simulations.

[0256] Since in practice no restraints are applied to the protein atoms in the DBC restraint, the correction term can be calculated in the gas phase. The free energy cost of adding the DBC to the ligand is computed through two steps: (a) the ligand is restrained by a flat-bottom restraint on its center of mass within a radius equal to DBCmax+1 A, and the free energy cost of this step is determined analytically; and (b) the ligand is restrained by the DBC restraint, with the free energy cost of this step calculated using thermodynamic integration (TI) via the Tinker and CoIvars libraries. The total correction from DBC is then obtained as the sum of these two contributions. To achieve better convergence, we ran the second step in two directions. In the first (referred to as forward), the ligand disappears over 202 values, while in the second (referred to as backward), the ligand reappears using the same 2 schedule. In theory, these two processes should converge to the same result. The final DBC correction is computed as the average of the forward and backward results.

[0257] In Table 2, we present the restraint correction terms derived from both DBC and Harmonic restraints. For DBC, the results are averaged over three replicas, with the errors corresponding to the standard deviation across the replicas.

[0258] Table 2: Amount of restraint from Distance-to-Bound-Conformation (DBC) and Harmonic Compound DBCcutoff (A) DBC (kcal / mol) Harmonic (kcal / mol)

[0259] 1 0.91 4.52 + 0.07 2.33

[0260] 2 0.79 5.00 + 0.05 2.51

[0261] 3 0.80 5.14 + 0.03 2.43

[0262] 4 0.91 3.77 + 0.05 2.33

[0263] 5 0.87 4.64 + 0.00 2.37

[0264] 6 0.86 5.47 + 0.02 2.38

[0265] 7 0.84 5.20 + 0.15 2.40

[0266] 8 0.88 4.88 + 0.02 2.36

[0267] 9 0.84 5.37 + 0.00 2.39

[0268] 10 0.96 4.19 + 0.04 2.28

[0269] 11 0.97 5.15 + 0.04 2.28

[0270] Lambda-ABF-OPES for Calculation of Absolute Binding Free Energies

[0271] Details of the Lambda-ABF and OPES-Explore methods can be found in the original publications (Invernizzi, M. et al., 2020; Invernizzi M. et al., 2022; Lagardere, L. et al., 2024) Here we provide a brief overview of their underlying theories.The Lambda-ABF algorithm adaptively computes the derivative of the free energy associated with the parameter using the thermodynamic integration (Tl) formula. This estimate is applied as a force directly to the simulated system, guiding the dynamics. As a result, the sampling of converges toward a uniform distribution, facilitating the crossing of free-energy barriers.

[0272] OPES-Explore, a collective-variable (CV)-based enhanced sampling technique, represents the latest advancement in the MetaDynamics family. It broadens the sampling of a system toward a target probability distribution, known as the well-tempered distribution. This approach employs Gaussian kernels to adaptively construct a bias potential, guiding the system’s exploration while preserving thermodynamic relevance. Critical parameters, such as the Gaussian kernel width and bias factor, are pivotal in defining the sampling’s scope and stability. Additionally, the barrier parameter (AE) ensures efficient transitions between basins while preventing access to irrelevant high-energy states.

[0273] In this implementation, the ABFE is calculated by continuously “alchemically” decoupling the ligand from its environment, both in complex with the protein and separately in bulk solvent. The standard free energy of binding is then determined using a thermodynamic cycle, which requires sampling of the alchemical Hamiltonians. We employed the Lambda-ABF approach as implemented in Tinker-HP / Colvars, in combination with OPES-Explore provided by the CoIvars library. This integration provides user-friendly convergence estimation without requiring postprocessing and allows seamless compatibility with other CV-based methods.

[0274] All simulations were performed at T = 300 K and P = 1 Atm using the BAOAB-RESPA integrator with a 3 fs time step under the NPT ensemble, employing four walkers to separately decouple the VDW and ELE interactions. More precisely, the polarizabilities and permanent multipoles of the ligand are scaled down to 0 in the electrostatic legs (ELE), and the van der Waals interactions between the atoms of the ligand and all the other ones are scaled down to 0 (leveraging softcore interactions) in the van der Waals legs (VD+W). In the complex phase, 30 ns were simulated for the VDW leg and 5 ns for the ELE leg. In the solvent phase, both the VDW and ELE legs were simulated for 5 ns each.

[0275] ELE and VDW decomposition of PMF:Here we present the electrostatic (ELE) and van der Waals (VDW) decomposition of the Potential of Mean Force (PMF). For each ligand, three replicas are reported.

[0276] Table 3: Energy component of PMF. ELE and VDW decomposition of the Potential of Mean Force (PMF) for the complex and solvent phases. For each ligand, three replicas are reported. The energies are presented in kcal / mol.

[0277] Complexation Solvation

[0278] Uul I IUUUI iu ncUUCd

[0279] H HELE VDW ELE VDW

[0280] 1 -173.44 61.22 -169.81 75.06

[0281] LI 2 -173.54 61.25 -169.48 75.18

[0282] 3 -173.35 61.28 -169.39 75.57

[0283] 1 -56.07 29.60 -56.26 45.91

[0284] L2 2 -56.13 29.40 -56.26 46.00

[0285] 3 -55.45 28.86 -56.16 45.54

[0286] 1 -80.42 22.81 -76.93 36.40

[0287] L3 2 -80.68 22.50 -76.88 36.19

[0288] 3 -80.49 22.51 -77.05 36.10

[0289] 1 -153.9 53.97 -151.62 68.24

[0290] L4 2 -154.78 55.11 -151.64 68.58

[0291] 3 -154.71 55.11 -151.31 68.50

[0292] 1 -36.49 10.32 -33.02 22.90

[0293] L5 2 -36.64 11.05 -32.83 23.10

[0294] 3 -36.40 11.55 -32.89 23.10

[0295] 1 -18.41 36.44 -17.35 50.25

[0296] L6 2 -18.54 36.42 -17.23 49.90

[0297] 3 -18.30 35.83 -17.14 50.02

[0298] 1 -68.73 34.38 -63.80 43.48

[0299] L7 2 -69.40 34.76 -63.83 43.17

[0300] 3 -69.30 33.69 -63.84 43.52

[0301] 1 -28.89 25.27 -27.61 36.79

[0302] L8 2 -28.58 23.96 -27.71 36.79

[0303] 3 -28.72 25.25 -27.74 37.09

[0304] 1 -28.30 8.06 -25.59 20.97

[0305] L9 2 -28.58 8.69 -25.61 21.37

[0306] 3 -28.96 8.84 -25.68 21.41

[0307] 1 -26.06 7.31 -22.53 17.68

[0308] L10 2 -27.42 9.13 -22.48 17.79

[0309] 3 -27.11 8.29 -22.51 17.97

[0310] 1 -7.44 4.76 -5.80 13.25

[0311] Lil, pose 1 2 -7.21 3.01 -5.81 13.01

[0312] 3 -7.41 2.93 -5.88 13.20

[0313] 1 -9.20 1.76 -5.77 13.25

[0314] Lil, pose 2 2 -9.16 1.51 -5.88 12.73

[0315] 3 -9.19 1.71 -5.85 12.75As mentioned earlier, for OPES-Explore, the barrier parameter should be sufficient to ensure that the method accurately captures the energy landscape and transition state barriers. We conducted tests using OPES-Explore alone with different barriers. The final AG values varied between runs, indicating that OPES-Explore alone struggles to compute the binding free energies of interest.

[0316] However, when combining OPES-Explore with Lambda-ABF, the system is primarily driven by Lambda-ABF, which facilitates the crossing of energy barriers, and also benefits from the rapid convergence of the Tl estimator. In this case, OPES-Explore mainly serves to push the system out of local minima. Therefore, the barrier parameter does not need to be equal to the free energy differences associated to the simulation. Our tests showed that setting the bias threshold of OPES-Explore to a minimum of 5 kcal / mol for both the ELE and VDW legs in the complex and solvent phases effectively accelerated the convergence across all ligands, regardless of their respective energy barriers. For ligand 11 , a threshold of 2 kcal / mol was sufficient for proper acceleration. Although a barrier of 5 kcal / mol also works well, since this ligand is a weak binder, the lower threshold of 2 kcal / mol ensures better stability in convergence. An adaptive sigma is used to determine the Gaussian kernel width, and the frequency for kernel deposition was set to 300 steps, which corresponds to 9 ps with a time step of 3 fs.

[0317] For Ligand 11 , during the 180 ns plain MD simulations, we observed two equally plausible binding modes in which the trifluorotoluene moiety was flipped by 90°, consistent with the result of prior art. Consequently, ABFE calculation was performed for each binding mode, resulting in a total of two binding free energy calculations. The results from multiple binding modes were combined into a single binding free energy value.

[0318] The free energy cost associated with the release of the DBC restraint is computed in the gas phase through Tl by progressively releasing it to a compatible harmonic distance restraint which can then be computed analytically.

[0319] Plain OPES setup for Asn140 dihedral FES

[0320] For the Apo state, we start the simulation from the X-ray structure (PDB ID: 4LYI). The same equilibration process as in the Holo state is applied here. Please refer to the main text for more details. After approximately 100 ns of production, we used thefinal frame for the plain OPES3 simulation. The biased simulation, using the dihedral as the collective variable (CV), is run for 160 ns with Tinker-HP and Plumed2. The BARRIER (AE) of OPES is set to 20 kJ / mol, with a PACE of 500 strides, which, with a time step of 10 fs, corresponds to 5 ns. The final reweighed FES is averaged over three blocks.

[0321] Well-Tempered Metadynamics (WTMD) and OPES-Explore

[0322] We ran a few tests to compare the performance of L-ABF-OPES with L-ABF-WTMD. For L-ABF-WTMD, the simulation is highly sensitive to the parameters of WTMD, such as the bias factor, hill height, and width of the Gaussian hills. Through trial and error, with the parameters set to a bias factor of 3000 K, hill weight of 0.1 , and hill width of 5.0, we were able to obtain reasonable results, with slightly faster convergence in L-ABF-OPES (see Fig. 6).

[0323] In the case of OPES-Explore alone, its exploratory nature, as discussed in Ref. 4, leads to inconsistencies across different repeats, especially when compared to L-ABF-OPES. Fig. 7 shows the FES of the VDW leg in the solvent phase of Ligand5 across three replicas. Each replica resulted in a different

[0324]

[0325] value.

[0326] In these simulations, we observed that the rotamers of Asn140 in the Apo state may adopt a conformation favorable in the Apo state but not in the Holo state (see Results and Discussion and SI for more details). This could lead to sampling the Holo state with an incorrect orientation of Asn140, potentially introducing artifacts. The restriction (through a restraint) to one of these rotamers as an Apo endpoint must be accounted for. Therefore, a positional restraint was applied to the heavy atoms of the Asn140 residue and its neighboring atoms (all heavy atoms of the backbone and those within 6 A of the ligand with a mild force constant of 2 kcal / mol / A2). This restraint was necessary to ensure that Asn140 and its neighboring atoms were properly sampled in both the Apo and Holo states, given the continuous switching of lambda between 0 (Apo) and 1 (Holo) in the Lambda-ABF-OPES simulations. Since both rotamers of Asn140 are equally favorable in the Apo state (see Fig. 3), an RTIn(2) correction was added to the final computed AG. In addition, to avoid artifacts during alchemical decoupling as described in SI, the Ccr atoms of BRD4 were restrained to their relaxed configuration.

[0327] Analytical Lambda Derivatives for Variational Many-Body PotentialsAs all Tl-based technique, Lambda-ABF requires the computation of potential derivatives with respect to the alchemical parameter A. We can use a simple interpolation of polarization between the end states which gives immediately the associated derivatives as the difference of these, but this is associated with an increase of the computational cost because of the need to solve 2 polarization equations at each timestep. Here we resorted to a more general formulation relying on the variational formulation of the many-body term and the Hellman-Feynman theorem: Epo / (r,A) = Epol(r, A, |i(r, A)) and

[0328]

[0329] because of the minimum conditions on the induced dipoles. Note that this formulation still holds for other many-body terms with similar variational formulation as is the case for fluctuating charges, continuum solvation models or QM-MM.

[0330] Results and Discussion

[0331] In the following section, we will first provide a brief overview of the general binding modes of the 11 ligands in the BRD4 system. We will then focus on the orientation of the Asn140 residue, a key residue that directly interacts with the ligands, in both the Apo and Holo states, highlighting any significant differences. Finally, we will examine the binding affinities and correlate these findings with the experimental values.

[0332] As depicted in Fig. 2(c), the 11 ligands are large, flexible, drug-like molecules, some of which are charged (Ligand 1 and 4). These characteristics make them an ideal set for evaluating the performance of the new method. All ligands target a common binding site. Asparagine (Asn140), located within the BC-loop binding pocket (Fig.

[0333] 2(a)), is the most critical residue, which interacts directly with the ligands. The side chain of Asn140 consists of an amide group (-CONH2) attached to a C / 3, which is itself connected to the Ccr of the backbone via rotatable bonds.

[0334] The 100 ns plain MD simulations for all ligands (except Ligand 11 , for which 180 ns was run) show that, in the Holo state, only one rotamer is favorable due to the formation of salt bridges between the ligands and the Asn140 residue (Fig. 2(b)). However, in the Apo state, both rotamer states can be present. To investigate this, we calculated the free energy associated with the two rotamers in the Apo state using the OPES method (see SI for more details). Fig. 3 shows that the two rotamers are equally favorable and that switching between them can happen naturally.In addition to Asn140, four conserved (polarizable) water molecules in the binding site of BRD4 also play a crucial role in stabilizing the ligands within the binding pocket, thereby enhancing binding affinities. As illustrated in Fig. 2(a), one of these water molecules forms a bridge between the Tyr97 residue and the ligand (except for Ligand 4). To preserve this stabilization, the X-ray water molecules were retained in the binding pocket during system preparation, as they contribute to a hydrogen-bond network involving the three other conserved water molecules and the protein.

[0335] Using the equilibrated structures, we carried out absolute binding free energy calculations employing the novel Lambda-ABF-OPES method. The correlation plot between experimental values and calculated results is shown in Fig. 4 and reported in Table 4. We observe a strong correlation between the experimental and calculated results, with a Pearson’s correlation coefficient (Pearson’ r) of 0.81. The analysis yielded a root mean square error (RMSE) of 1.1 kcal / mol and a mean absolute error (MAE) of 0.9 kcal / mol. These results are consistent with experimental data and align well with those reported in prior art. The key advantage of this method lies in its user-friendliness, reduced computational resource requirements, and rapid convergence. Compared to Lambda-ABF, the traditional fixed Lambda approach, and the newly developed CV-based approach, this method offers significantly improved computational efficiency while maintaining high accuracy. These attributes make it particularly well-suited for large-scale or high-throughput applications.

[0336] To evaluate the robustness of the method, we performed a detailed convergence test by analyzing the convergence time of AG and comparing the results with those obtained using Lambda-ABF alone. The convergence time was defined as the point from which all subsequent data points remain within the specified tolerance of the mean, providing a clear metric for determining when the simulation data stabilizes. A tolerance of 0.2 kcal / mol was used in this analysis, representing just 20% of the commonly accepted convergence threshold for binding free energy calculations. This strict criterion underscores the precision and reliability of the method. Figure 5 illustrates an example of a AG convergence plot over time for the ELE and VDW legs of the Ligand 8 complex phase using Lambda-ABF and Lambda-ABF-OPES method in a single replica. In this case, the acceleration in convergence time is a factor of 9for ELE and 4 for VDW. When considering different replicas, the convergence speedup for ELE ranges from 5 to 9 times, while for VDW, it ranges from 3 to 5 times.

[0337] The convergence analysis results for all ligands using Lambda-ABF-OPES are summarized in Table 4. The average convergence time for the complex phase was less than 3 ns for the ELE leg and approximately 20 ns for the VDW leg, highlighting the rapid stabilization of the AG in these components. In the solvent phase, the ELE leg showed convergence for most ligands within 1-2 ns or less, while the VDW leg converged within 2-3 ns. These results demonstrate a marked improvement in convergence speed compared to other currently used methods, which typically

[0338] requires significantly longer simulation times to achieve comparable stability.

[0339] In practical applications, such efficiency gains can translate into a substantial reduction in computational costs, enabling more extensive exploration of chemical space or greater statistical sampling within the same resource constraints. The combination of simplicity, accuracy, and computational efficiency makes this approach a promising tool for the drug discovery field.

[0340] Table 4: Summary of the BRD4 binding free energy results using Lambda-ABF-OPES. The experimental (AGexp) and calculated (AGcalc) values for each ligand are presented. All AG values are reported in kcal / mol. The calculated AGcalc and associated errors represent the mean and standard error of the mean, derived from three replicates for each ligand. The PDB files used as input are listed. The average convergence time over three replicas for each ELE and VDW leg in the complex and solvent phases is reported in ns per replica. The total simulation time for the ELE and VDW legs in the complex phase is 5 ns and 30 ns per walker, respectively. For the solvent phase, each leg is run for 5 ns per walker.

[0341] ns per walker / \ fir fli "

[0342] Compound AGexp AGcalc AGexp PDB Complex Solvent ELE VDW ELE VDW LI -9.8 ± 0.1 -10.56 + 0.46 -0.85 4OGI 0.9 22.5 1.9 3.5 L2 -9.6 ±0.1 -8.26 + 0.22 1.34 3MXF 2.2 19.2 0.9 2.2 L3 -9.0 + 0.1 -9.19 + 0.30 -0.19 4MR3 2.6 17.8 0.4 2.4 L4 -8.9 + 0.1 -10.13 + 0.15 -1.23 4OGJ 2.6 24.8 2.7 3.5 L5 -8.8 + 0.1 -8.27 + 0.52 0.09 4J0R 2.1 15.9 0.9 2.7 L6 -8.2 + 0.1 -6.75 + 0.30 1.45 3U5L 0.6 20.2 0.5 3.5 L7 -7.8 + 0.1 -6.42 + 0.65 1.38 4MR4 0.9 20.0 1.4 3.4 L8 -7.4 + 0.1 -5.46 + 0.53 1.94 3U5J 0.7 20.7 0.4 2.5 L9 -7.3 + 0.1 -7.57 + 0.12 -0.27 3SVG 1.5 17.4 0.5 2.7LIO -6.3 ± 0.0 -7.05 + 0.33 -0.75 4HBV 1.8 14.8 0.2 1.6 Lil -5.6 -4.98 + 0.58 0.62 Model 0.2 10.7 0.2 1.9

[0343] Conclusions

[0344] Accurately predicting the binding affinity between small molecules and their target proteins remains a critical challenge in drug discovery, with far-reaching implications for the speed and efficiency of therapeutic development.

[0345] Here, we introduced a new hybrid approach that combines Lambda-ABF with the exploratory version of the OPES method. This novel integration leverages the complementary strengths of ABF and OPES to overcome critical limitations in alchemical free energy calculations, including inefficient exploration of configurational space and kinetic trapping in energy landscapes. By applying biases to the Lambda collective variable (CV) and incorporating the AMOEBA polarizable force field alongside DBC restraints, our method achieves unprecedented levels of sampling efficiency, up to nine times faster than the original Lambda-ABF technique.

[0346] Our application of this hybrid method to a diverse set of 11 drug-like molecules targeting BRD4 bromodomains yielded close alignment between our computational results and experimental data with a mean absolute error of 0.9 kcal / mol. Importantly, this was achieved while using the highly accurate AMOEBA polarizable force field, demonstrating the feasibility of this approach for real-world drug discovery applications. This approach can be naturally extended to neural networks methodologies including Machine Learning Interatomic Potentials (MLIP) and Foundation models whose additional computational cost for free energy computations compared to FFs has limited to date their use in production.

[0347] By integrating state-of-the-art methodologies and harnessing their synergistic advantages, this work provides a robust tool for the rapid and reliable advancement of novel therapeutics. Though applied in an alchemical context in this implementation, this methodology shows promise for broader applicability within enhanced sampling techniques and lays the groundwork for its integration into a more general framework, which we plan to further investigate in future work.

[0348] References

[0349] Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo freeenergy estimation: Umbrella sampling. J. Comput. Phys. 1977, 23, 187-199.

[0350] Darve, E.; Pohorille, A. Calculating free energies using average force. J. Chem. Phys. 2001, 115, 9169-9183.

[0351] Darve, E.; Rodnguez-Gomez, D.; Pohorille, A. Adaptive biasing force method for scalar and vector free energy calculations. J. Chem. Phys. 2008, 128.Comer, J.; Gumbart, J. C.; Henin, J.; Lelievre, T.; Pohorille, A.; Chipot, C. The adaptive biasing force method: Everything you always wanted to know but were afraid to ask. J. Phys. Chem. B 2015, 119, 1129-1151.

[0352] Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. 2002, 99, 12562-12566.

[0353] Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 826-843.

[0354] Invernizzi, M.; Parrinello, M. Rethinking Metadynamics: From bias potentials to probability Distributions. J. Phys. Chem. Lett. 2020, 11, 2731-2736.

[0355] Invernizzi, M.; Parrinello, M. Exploration vs convergence speed in adaptive-bias enhanced sampling. J. Chem. Theory Comput. 2022, 18, 3988-3996.

[0356] Maragliano, L.; Vanden-Eijnden, E. A temperature accelerated method for sampling free energy and determining reaction pathways in rare events simulations. Chem. Phys. Lett. 2006, 426, 168-175.

[0357] Henin, J.; Lelievre, T.; Shirts, M. R.; Valsson, O.; Delemotte, L. Enhanced Sampling Methods for Molecular Dynamics Simulations [Article v1.0]. Living J. Comput. Mol. Sci. 2022, 4, 1583. Lagardere, L.; Maurin, L.; Adjoua, O.; El Hage, K.; Monmarche, P.; Piquemal, J.- P.; Henin, J. Lambda-ABF: Simplified, portable, accurate, and cost-effective alchemical free-energy computation. J. Chem. Theory Comput. 2024, 20, 4481-4498, PMID: 38805379.

[0358] Salari, R.; Joseph, T.; Lohia, R.; Henin, J.; Brannigan, G. A Streamlined, general approach for computing ligand binding free energies and its application to GPCRbound cholesterol. J. Chem. Theory Comput. 2018, 14, 6560-6573, PMID: 30358394.

[0359] Fiorin, G.; Klein, M. L.; Henin, J. Using collective variables to drive molecular dynamics simulations. Mol. Phys. 2013, 111, 3345-3362.

[0360] Fiorin, G.; Marinelli, F.; Forrest, L. R.; Chen, H.; Chipot, C.; Kohlmeyer, A.; Santuz, H.;

[0361] Henin, J. Expanded functionality and portability for the CoIvars library. J. Phys. Chem. B 2024, 128, 11108-11123.

Claims

Claims1. A computer-implemented method (100) for computing aft free-energy difference between two end-state Hamiltonians of a molecular system by molecular dynamics simulation, the method (100) comprising:- a) obtaining (110) at least one collective variable describing the molecular system; - b) During the molecular dynamics simulation, applying a biasing scheme (120) to the at least one collective variable, said a biasing scheme being updated on the fly, said biasing scheme comprising:o an adaptive biasing force (121) configured to drive sampling of the least one collective variable toward a flattened distribution, ando a bias potential (122) added as an additional potential energy term, the bias potential being defined as a function of a target probability distribution of the least one collective variable, wherein the target probability distribution is a well-tempered distribution, and wherein the bias potential is constructed using Gaussian kernels; and- c) computing the free-energy difference (130) by thermodynamic integration of energy derivatives calculated from the collected information, with respect to at least one collective variable, obtained during the step b).

2. The computer-implemented method (100) according to claim 1 , wherein at least one collective variable is selected among alchemical variable, distance, angle, dihedral, root mean square deviation, coordination number and / or machine-learning based collective variables.

3. The computer-implemented method (100) according to claim 1 or 2, wherein the at least one collective variable comprise an alchemical parameter A that interpolates between the two end-state Hamiltonians and the biasing scheme is applied to the alchemical parameter A; the adaptive biasing force being configured to drive sampling of the alchemical parameter A toward a flattened distribution and the bias potential being added as an additional potential energy term, the bias potential being defined as a function of a target probability distribution of the alchemical parameter A.

4. The computer-implemented method (100) according to anyone of claims 1 to 3, wherein, at the step b), a weighted kernel density estimation is employed.

5. The computer-implemented method (100) according to claim 4, wherein the weighted kernel density estimation uses an adaptive kernel width that decreases as the number of effective samples increases.

6. The computer-implemented method (100) according to any one of claims 1 to 5, wherein, at the step b), it comprises a step of estimating an underlying probability distribution of the CV and a step of adjusting the bias potential so that the simulated distribution of the CVs approaches the target probability distribution.

7. The computer-implemented method (100) according to claim 6, wherein estimating the underlying probability distribution comprises performing a weighted kernel density estimation of the CV distribution using simulation samples, weighting each sample by a factor related to the current bias to approximate the unbiased probability; preferably it also comprises compressing the kernel representation by merging kernels in previously sampled regions to limit computational growth.

8. The computer-implemented method (100) according to anyone of claims 1 to 7, wherein, at the step b), the adaptive biasing force is computed from an on-the-fly estimate of the equilibrium (unbiased) target probability distribution of the molecular system.

9. The computer-implemented method (100) according to anyone of claims 1 to 8, wherein, at the step b), the adaptive biasing force is computed from a biasing potential issued from an on-the-fly estimate of the current biased probability distribution of the molecular system.

10. The computer-implemented method (100) according to anyone of claims 1 to 9, wherein, at the step b), the adaptive biasing force is computed from a biasing potential issued from an on-the-fly estimate of the well-tempered distribution of the molecular system being sampled.

11. The computer-implemented method (100) according to anyone of claims 1 to 10, wherein it is a computer-implemented method for computing alchemical binding free- energy in the molecular system.44 / 5112. The computer-implemented method (100) according to anyone of claims 1 to 11 , wherein the at least one collective variable describing the molecular system is an alchemical parameter A which is between 0 and 1.

13. The computer-implemented method (100) according to claim 12, wherein it comprises a step of defining reflective boundary conditions for A=0 and A=1 , ensuring that said alchemical parameter remains within [0,1] during the simulation.

14. The computer-implemented method (100) according to anyone of claims 1 to 13, wherein the Gaussian kernels are defined with parameters comprising at least one from: Gaussian kernel width and barrier parameter; preferably barrier parameter.

15. The computer-implemented method (100) according to claim 14, wherein the barrier parameter (AE) is selected to allow efficient transitions between basins while preventing access to irrelevant high-energy states.

16. The computer-implemented method (100) according to anyone of claims 14 or 15, wherein an adaptive sigma is used to determine the Gaussian kernel width.

17. The computer-implemented method (100) according to anyone of claims 14 to 16, wherein a frequency for kernel deposition is set as at least one of: a number of simulation steps and a deposition time interval.

18. The computer-implemented method (100) according to anyone of claims 1 to 17, wherein the step b) is applied concurrently by multiple-walker sampling instances of the at least one collective variable describing the molecular system.

19. The computer-implemented method (100) of the previous claim, wherein the multiplewalker sampling instances exchange or share adaptive biasing force data so as to collectively explore different regions of phase space.

20. The computer-implemented method (100) of claim 18, wherein each multiple-walker sampling instances periodically transmits locally accumulated biasing data to a central data structure or to all other multiple-walker sampling instances.

21. The computer-implemented method (100) of claim 18, the at least one collective variable describing the molecular system is an alchemical parameter A and wherein each multiple-walker sampling instances :updates the molecular system’s atomic coordinates and velocities, updates the alchemical lambda parameter A,- calculates a local mean force on alchemical lambda parameter A and adjusts the adaptive bias, andperiodically exchanges mean force or bias information with other multiplewalker sampling instances;o aggregating the simulation data from all multiple-walker sampling instances to compute a free-energy difference for each value of A; and o integrating the aggregated mean force data to obtain the alchemical free-energy of the molecular system.

22. The computer-implemented method (100) of claim 18, the at least one collective variable describing the molecular system is an alchemical parameter A and wherein each multiple-walker sampling instances :updates the molecular system’s atomic coordinates and velocities, updates the alchemical lambda parameter A,- calculates a local mean force on alchemical lambda parameter A and adjusts the adaptive bias including the use of a bias potential expressed as a function of a target probability distribution of the molecular system, and periodically exchanges mean force or bias information with other multiplewalker sampling instances;o aggregating the simulation data from all multiple-walker sampling instances to compute a free-energy gradient for each value of A; and o integrating the aggregated mean force data to obtain the alchemical free-energy of the molecular system.

23. The computer-implemented method (100) according to anyone of claims 1 to 22, wherein the molecular system comprises a ligand, the at least one collective variable describing the molecular system is an alchemical parameter A and A interpolates between a fully coupled state and a decoupled state of the ligand.

24. The computer-implemented method (100) according to claim 23, wherein a binding free energy is calculated by separately decoupling van der Waals interactions and electrostatic interactions.

25. The computer-implemented method (100) according to claim 24, wherein polarizabilities and permanent multipoles of the ligand are scaled down to 0 in electrostatic legs, and van der Waals interactions between atoms of the ligand and all other atoms are scaled down to 0 in van der Waals legs.

26. The computer-implemented method (100) according to claim 25, wherein the van der Waals legs leverage softcore interactions.

27. The computer-implemented method (100) of anyone of claims 1 to 26, wherein it further comprises leveraging Distance-to-Bound-Configuration (DBC) restraints to at least a part of the molecular system; preferably a ligand is partially restrained by a distance-to-bound-configuration (DBC) variable during equilibration.

28. The computer-implemented method (100) according to claim 27, wherein the distance-to-bound-configuration (DBC) variable is defined as a root mean square deviation of ligand atoms for each frame, with an alignment performed relative to atoms of a receptor binding site.

29. The computer-implemented method (100) according to claim 28, wherein a DBC cutoff is selected as a DBC value within a 95% interval of a DBC distribution monitored during plain molecular dynamics of the ligand, and the selected DBC cutoff is used as a cutoff for the DBC restraint.

30. The computer-implemented method (100) according to claim 29, wherein a flat- bottomed harmonic restraint is applied to the DBC above the DBC cutoff.

31. The computer-implemented method (100) according to claim 30, wherein a free energy cost associated with a release of the DBC restraint is computed in a gas phase through thermodynamic integration by progressively releasing the DBC restraint to a compatible harmonic distance restraint which is computed analytically.

32. The computer-implemented method (100) according to any of claims 27 to 31 , wherein selecting ligand atoms used for the DBC variable comprises monitoring a root mean square fluctuation of ligand heavy atoms over a plain molecular dynamics simulation and selecting atoms based on the root mean square fluctuation.

33. The computer-implemented method (100) according to anyone of claims 1 to 32, the at least one collective variable describing the molecular system is an alchemical parameter A and wherein A is a coupling parameter controlling the interpolation between two end-state Hamiltonians (e.g., a ligand fully non-interacting vs fully interacting).

34. The computer-implemented method (100) of anyone of claims 1 to 33, the at least one collective variable describing the molecular system is an alchemical parameter A and wherein performing thermodynamic integration comprises continuously integrating an ensemble-averaged derivative of the molecular system’s potential energy with respect to the alchemical parameter A over a range of A values to obtain the alchemical free-energy of the molecular system.

35. The computer-implemented method (100) of anyone of claims 1 to 34, wherein the molecular simulation employs a force field selected from: fixed-charge force fields (e.g. AMBER, CHARMM, OPLS-AA), machine learning force fields, coarse-grained force fields, reactive force fields, bond-order based force fields, ionic force fields and / or polarizable force fields (e.g. AMOEBA); or the molecular simulation employs an Hamiltonian obtained through Quantum Mechanical / Molecular Mechanical (QM / MM) methods, and the lambda-dynamics adaptive biasing force approach is applied under either one or several types of force-field models.

36. The computer-implemented method (100) of anyone of claims 1 to 35, the at least one collective variable describing the molecular system is an alchemical parameter A and wherein the thermodynamic integration comprises numerically integrating mean forces exerted along said alchemical parameter A, thereby obtaining a free-energy difference between a bound state and an unbound state.

37. The computer-implemented method (100) according to anyone of claims 1 to 36, the at least one collective variable describing the molecular system is an alchemical parameter A and wherein free-energy is monitored in real time by evaluating partial integrals of a mean force with respect to A and analyzing resulting estimates over cumulative simulation time.

38. The computer-implemented method (100) of anyone of claims 1 to 37, wherein it comprises the computation of potential derivatives with respect to the alchemical lambda parameter through interpolation of polarization between the end states which48 / 51gives immediately the associated derivatives as the difference of these (for the polarizable force field).

39. The computer-implemented method (100) of anyone of claims 1 to 38, the at least one collective variable describing the molecular system is an alchemical parameter A and wherein step (c) comprises computing a free-energy profile as a function of A by integrating said mean forces with respect to A, and obtaining the alchemical free- energy difference between the two end-state Hamiltonians from said integrated profile.

40. The computer-implemented method (100) of anyone of claims 1 to 38, wherein step (c) is carried out during step (b) so as to provide an on-the-fly estimate of the alchemical free-energy difference during the simulation.

41. The computer-implemented method (100) of anyone of claims 1 to 38, the at least one collective variable describing the molecular system is an alchemical parameter A and wherein step (c) comprises monitoring free-energy in real time by evaluating partial integrals of the mean force with respect to A and analyzing resulting estimates over cumulative simulation time.

42. The computer-implemented method (100) of anyone of claims 1 to 38, the at least one collective variable describing the molecular system is an alchemical parameter A and wherein step (c) comprises integrating mean forces over a range of A values sampled continuously in a single A-dynamics simulation, without breaking the A range into discrete fixed-A windows.

43. The computer-implemented method (100) of anyone of claims 1 to 38, the at least one collective variable describing the molecular system is an alchemical parameter A and wherein step (c) comprises aggregating mean-force data obtained from a plurality of multiple-walker sampling instances and integrating the aggregated meanforce data with respect to A to obtain the alchemical free energy.

44. The computer-implemented method (100) of claim 43, the at least one collective variable describing the molecular system is an alchemical parameter A and wherein aggregating comprises computing a free-energy gradient (mean-force profile) for each value of A from the plurality of multiple-walker sampling instances, and integrating said free-energy gradient with respect to A.49 / 5145. The computer-implemented method (100) of anyone of claims 1 to 38, wherein a standard free energy of binding is determined using a thermodynamic cycle based on (i) a complex phase in which the ligand is decoupled in complex with a protein and (ii) a solvent phase in which the ligand is decoupled in bulk solvent, each phase’s decoupling free energy being computed by step (c).

46. The computer-implemented method (100) of anyone of claims 1 to 38, the at least one collective variable describing the molecular system is an alchemical parameter A and wherein the mean forces integrated in step (c) comprise potential-energy derivatives with respect to A computed during the simulation.

47. The computer-implemented method (100) of claim 46, wherein the molecular simulation employs a polarizable force field, and wherein the potential-energy derivatives comprise a derivative contribution associated with polarization energy.

48. The computer-implemented method (100) of anyone of claims 1 to 47, wherein it is used to: compute a binding free energy of a ligand to a target; compute an absolute binding free energy of a ligand; compute a binding affinity between a ligand and a receptor; compute a free-energy profile for protein folding; compute free-energy differences between molecular conformational states; compute activation free energies; compute reaction free energies; compute solvation free energies; compute partition coefficients; compute dissociation constants; compute association constants; compute relative stabilities of molecular conformers; compute free-energy differences in enzyme-substrate complexes; compute free energies for conformational transitions; compute free energies for phase transitions; compute free energies of adsorption; compute free energies for ligand-induced conformational changes; compute free energies in molecular docking simulations; compute free energies for applications in drug design; compute free energies for stability prediction of molecular systems; compute free energies in polymer systems; compute free energies in solid- state systems; compute free energies in solvated systems; compute free energies for protein-ligand interactions; compute free energies for nucleic acid-ligand interactions; compute free energies in molecular dynamics simulations; compute free energies for conformational sampling; and / or compute free energies for estimating reaction rates.50 / 5149. A system for computing free-energy surface in a molecular system, the system comprising one or several processors and one or several non-transitory computer- readable memories storing instructions which, when executed by the processor(s), cause the processor(s) to perform the method according to claims 1 to 47.

50. A non-transitory computer-readable medium storing instructions which, when executed by a processor, cause the processor to perform a computer-implemented method for computing a free-energy surface of a molecular system according to anyone of claims 1 to 47.

Images