Atomistic Model for Simulations of the Sedative Hypnotic Drug 2,2,2-Trichloroethanol.

2,2,2-Trichloroethanol (TCE) is the active form of the sedative hypnotic drug chloral hydrate, one of the oldest sleep medications in the market. Understanding of TCE's action mechanisms to its many targets, particularly within the ion channel family, could benefit from the state-of-the-art computational molecular studies. In this direction, we employed de novo modeling aided by the force field toolkit to develop CHARMM36-compatible TCE parameters. The classical potential energy function was calibrated targeting molecular conformations, local interactions with water molecules, and liquid bulk properties. Reference data comes from both tabulated thermodynamic properties and ab initio calculations at the MP2 level. TCE solvation free energy calculations in water and oil reproduce a lipophilic, yet nonhydrophobic, behavior. Indeed, the potential mean force profile for TCE partition through the phospholipid bilayer reveals the sedative's preference for the interfacial region. The calculated partition coefficient also matches experimental measures. Further validation of the proposed parameters is supported by the model's ability to recapitulate quenching experiments demonstrating TCE binding to bovine serum albumin.

Introduction

Chloral hydrate is a sedative hypnotic drug used for short-term insomnia treatment. First synthesized in 1832,[1] chloral hydrate is the oldest sleep medication in the market. Chloral hydrate has been widely used as sedative in children undergoing clinical procedures. Its prescription in pediatrics is recommended for certain diagnostic procedures such as neurological imaging,[2−5] echocardiography,[6] and auditory brainstem response testing[7] when the patients do not respond to other agents.[8]

After oral administration, the chloral hydrate prodrug is rapidly converted into 2,2,2-trichloroethanol (C2H3Cl3O2). This active metabolite acts at the barbiturate recognition site on γ-aminobutyric acid (GABAA) receptors, eliciting sedation.[9] Whole-cell patch-clamp recording has shown that 2,2,2-trichloroethanol (TCE) potentiates GABAA-activated chloride current in mouse hippocampal neurons,[10] enhancing synaptic transmission. In addition to GABAA receptors, electrophysiology essays reported that TCE can interact with several other molecular targets. For example, TCE inhibits excitatory N-methyl-d-aspartate receptor and kainate-activated currents in mouse hippocampal neurons at clinical concentrations, more potently than ethanol.[11] Recently, it has also been shown that TCE modulates the recombinant human two-pore-domain potassium channels TREK-1 and TRAAK in a reversible, concentration-dependent manner.[12] Furthermore, tryptophan fluorescence quenching experiments demonstrated that TCE binds to bovine (BSA) and human (HSA) serum albumin plasma proteins.[13]

Despite the fact that TCE potentiates a range of protein receptors, its action mechanisms at the microscopic level are not clear. In the context of molecular dynamics (MD) simulations, understanding how TCE modulates these proteins requires an accurate atomistic model for both receptor and ligand. Although high-resolution crystallographic structures for some of these targets are available in the Protein Data Bank, an all-atom TCE model is still missing. We therefore present TCE parameters compatible with the CHARMM additive force field for biomolecules, useful for MD simulations.[14] The model is based on target quantities including molecular conformations, bulk phase properties, partition coefficient in water and oil, and TCE’s binding affinity to BSA.

Results and Discussion

TCE Parameters

TCE presents two main conformers, trans and gauche (Figure ), with electronic properties calculated at the MP2 level shown in Table . The magnitude of the trans conformation dipole moment |μ⃗| differs significantly from that of the gauche one (3.2 D against 1.6 D, respectively). Despite that, the isotropic portion of the polarizability tensor α̅ (au) shows almost the same magnitude for both geometries (53.0 and 52.1 au, respectively), with all individual components being very similar. The values for nondiagonal components suggest that both TCE conformers display an isotropic polarizable character. Gauche is the most stable at the MP2 level; the energy difference between geometries is 3.22 kcal/mol. All-atom parameters were thus calibrated by taking the gauche conformer as the molecular target.

Figure 1

Ball-and-stick representation of TCE conformers (A) gauche- and (B) trans-optimized at the MP2 level. Atom names are indicated. Graphics were rendered using GaussView.[15]

Table 1

Electronic Properties of trans and gauche Conformers of TCE Calculated at the MP2 Level

	|μ⃗| (D)	α̅ (au)	axx	axy	ayy	axz	ayz	azz
trans	3.2	52.1	53.6	0.0	53.6	–4.2	0.0	49.1
gauche	1.6	53.0	55.5	–0.7	55.6	–5.9	0.3	47.7

The dipole moment of TCE at the quantum mechanical (QM) level is overestimated by 43% in the mechanical molecular (MM) model, preserving vector orientation.[14] MM parameters for bond and valence angles reproduce QM geometry within an acceptable margin of error with deviations up to 0.03 Å and 3°, respectively.[14] TCE intramolecular interactions also agree with same data obtained from large-angle X-ray scattering experiments.[16] Torsional angles Cl–C2–C1–O and H–O–C1–C2 are used as conformational descriptors, the latter distinguishing trans and gauche conformers. As shown in Figure , their MM potential energy surfaces (PESs) fit QM data within 0.125 kcal/mol, comparable to kBT. Parameters and associated errors for partial charges and Lennard-Jones (LJ) and bond parameters are presented as Supporting Information Tables S1–S3. Note that Lennard-Jones parameters used to treat van der Waals interactions of TCE were assigned by transferability of similar chemical types available at the CHARMM force field (cf. Computational Methods).

Figure 2

Potential energy surfaces for torsional angles Cl–C2–C1–O (red) and H–O–C1–C2 (blue). Data is obtained from fully relaxed torsional scans at the MP2 level (solid lines) and empirical model (dashed lines). Note the two local minima for H–O–C1–C2 torsional angle at −70 and +70°, both gauche conformers.

Pure Solvent Properties

As shown in Table , calculated density and enthalpy of vaporization of TCE agree with experimental reference values, reproducing physicochemical properties of the bulk phase. Radial pair distribution functions (RDFs) as computed from pure solvent MD simulations of TCE show well-resolved peaks at 1.79, 2.72, 2.92, 3.12, and 3.94 Å, respectively, assigned to bond (−) and nonbonded (···) interactions C2–Cl1, C1···Cl1, Cl1···Cl3, O1···Cl1, and O1···Cl3 (Figure ). A hydrogen bond between O1 and O1 is also resolved at 2.92 Å. All peaks coincide with the reference RDF description for pure TCE, indicating that intra- and intermolecular interactions are correctly reproduced at the atomic level.[16] Note that simulation of TCE is predominantly populated by the gauche conformer with an occupancy probability of 78%, as expected from its larger stability. The gauche conformer thus accounts for most of the nonbonded interactions resolved in RDF analysis.

Table 2

Bulk Phase Properties of TCE

	calculated	reference	deviation (%)
density (g/mL)	1.59 ± 0.013	1.49a	6.0
enthalpy of vaporization (kcal/mol)	11.42 ± 0.280	10.82b	5.5

Experimental value available at ACS.

Theoretical value available at ACS. Calculated using Advanced Chemistry Development (ACD/Labs) Software V11.02 (1994-2018 ACD/Labs).

Figure 3

Solvation properties of TCE. First and second columns show respectively molecular systems and corresponding radial pair distribution functions (RDFs) for (a, b) pure TCE used in MD bulk phase simulations; (c, d) TCE in water; and (e, f) TCE in hexane. RDF of TCE’s oxygen O1 relative to water hydrogen atoms evidences a first solvation layer at approximately 1.8 Å (d), whereas RDF of TCE in hexane shows a first solvation layer between 5 and 6 Å (f). (g) Licorice representation of the TCE gauche conformer. Free energy change as function of TCE’s coupling/decoupling parameter λ in (h) bulk water phase and (i) in hexane. Error bars were estimated by the simple overlap sampling (SOS) method.[18] Molecular images were rendered using visual MD (VMD).[19]

Solvation Properties

TCE’s solvation free energy into hexane (oil) was calculated via free energy perturbation (FEP) and amounts to ΔGoil = −4.60 ± 0.001 kcal/mol, in agreement with TCE’s expected apolar character.[17] In spite of its favorable interaction with the apolar medium, TCE’s solvation free energy in water is within the same range as that in oil (ΔGwater = −4.83 ± 0.001 kcal/mol), suggesting that TCE is a lipophilic yet nonhydrophobic molecule. Figure displays RDF profiles and solvation free energy for TCE in both water and oil (hexane).

Transfer of TCE Across Membrane

To characterize TCE’s partition into the membrane, adaptive biasing force (ABF) simulations were employed to compute the potential of mean force (PMF) of a single TCE crossing a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer in the transmembrane z direction (Figure ). Analysis of PMF evidences a clear energy minimum at the phospholipids’ headgroup–tail interface region, with a stabilization of ΔGPOPC° = −4.32 ± 0.690 kcal/mol relative to bulk. TCE behavior accompanies the trend of some anesthetics such as isoflurane and sevoflurane, which show a distinct preference for the membrane interface.[20] For comparison, Figure also shows sevoflurane PMF, calculated under the same conditions as for TCE (anesthetic parameters come from Barber et al.[21]). PMF-derived atomic density for TCE is depicted in Figure c.

Figure 4

Transfer of TCE across membrane. (a) TCE–membrane system considered in ABF simulations. Phospholipid tails (i), membrane headgroups (ii), and water (iii) are indicated. (b) Potential of mean force (PMF) of TCE along the transmembrane direction. Accumulated statistical error for the PMF was estimated following eq (cf. Computational Methods) and amounts to 0.69 kcal/mol, indicating convergence of calculations. Sevoflurane PMF is shown for comparison. (c) Atomic density of TCE as computed from PMF in (b). (d) Atomic density profiles for water and membrane moieties. (e) Pair correlation between TCE and lipid headgroups. The first peak at 1.5 Å configures a typical hydrogen bond. The inset shows the molecular view of TCE binding to membrane headgroups via hydrogen bonds. Molecular images were rendered using VMD.[19]

It is worth noting that TCE is a halogenated molecule containing a hydroxyl group that can act as a hydrogen bond donor. Indeed, TCE–membrane headgroup RDF indicates a first peak at about 1.5 Å, configuring a typical hydrogen bond between the ligand and the lipid (Figure e). The hydroxyl group observed in TCE is absent in the general anesthetic sevoflurane, thus rationalizing their differential stabilization within the membrane core. Although the former displays modest stabilization (less than 1 kcal/mol) in water, the energy difference in the latter favors membrane core by approximately 3 kcal/mol.

The PMF of TCE was integrated to estimate its POPC–water partition coefficient (log K), eq . The calculated and tabulated theoreticala values, respectively, log K = 1.69 ± 0.001 and 0.97 ± 0.389 at 298.15 K, are within the same order of magnitude, supporting model’s ability to reproduce ligand’s lipid solubility properties.

Dissociation Constants for BSA Sites

TCE is known to interact and bind to both human and bovine serum albumin (HSA and BSA, respectively). To further validate the proposed TCE parameters, model’s affinity to these binding partners are confronted with ones determined from fluorescence quenching experiments.[13]

Briefly, in the experiment, BSA tryptophan residues W134 and W213 were effectively quenched by TCE at pH 7.0, allowing to pin down an apparent TCE dissociation constant of KD = 3.3 ± 0.3 mmol/L. HSA contains a single tryptophan residue, namely, W214, which is analogous to BSA’s W213. Thus, a combination of HSA and BSA quenching experiments meant individual dissociation constants of TCE from sites in the vicinity of BSA’s W134 and W213 could be estimated, respectively, as KD = 2.1 ± 0.1 and 12.0 mmol/L.[13]

Because quenching of W134 and W213 suggests two BSA binding sites for TCE, here labeled as s1 and s2, molecular docking and free energy perturbation (FEP) calculations were performed with the aim at reproducing experimental data (cf. Computational Methods) (Figure ). All measures (calculated and experimental) are reported in Table . The progression of free energy changes as a function of λ parameter for BSA binding sites during the course of five successive coupling/decoupling simulations is shown in Figure . Curves are well-behaved, which indicates that the time scale was small enough to correlate multiple homologous windows. As shown in Table , our predicted dissociation constants for ligand binding at sites s1 and s2 are in the same order of magnitude as that of the experimental estimates. Besides, the predicted constants also agree qualitatively with measurements as the experimentally determined higher affinity of TCE to site s1 is recapitulated in the calculations.

Figure 5

TCE binding sites to BSA. (a) Ensemble-average structure of BSA (cartoon), along with tryptophan residues (yellow licorice) and set of TCE centroid configurations (orange points) determined from docking searches. (b, c) Atomistic details of TCE within binding sites s2 and s1, respectively, composed of residues K20, V40, K131, W134, G135 and L197, W213, S343, L346, S453, L480, V481. Site residues are represented by surfaces and colored by physical and chemical properties. Molecular images were rendered using VMD.[19]

Table 3

Free Energy Changes for BSA Binding Sites s1 and s2 and Their Respective Dissociation Constantsa

	k	W*(SOS)	μ̅(SOS)	ΔGBSA°	KD (calc)	KD (exp)
s1	1.078	–11.60 ± 0.003	–4.83 ± 0.004	–6.77 ± 0.007	2.79c	2.1 ± 0.1e
s2	0.013	–7.60 ± 0.007	–4.83 ± 0.004	–2.77 ± 0.010	35.82c	12.0e
s1, s2				–5.78b	99.94d

Units for k, W*, μ̅, ΔGBSA°, and KD are kcal/mol/Å2, kcal/mol, kcal/mol, kcal/mol, and mmol/L, respectively.

Calculated standard free energy of binding a single ligand to both binding sites (independent events).

Calculated dissociation constant values.

Calculated aggregate dissociation constant value.

Experimental dissociation constant values from ref (13).

Figure 6

Computed free energy changes W* as a function of decoupling parameter λ for BSA binding sites. Error bars were calculated using the simple overlap sampling (SOS) method.[18]

Sites s1 and s2 are both buried in the protein and partially accessible to solvent. Independent MD simulations of the protein with bound ligands at sites s1 and s2 show that TCE remains confined to the binding sites, hydrated by few water molecules and in close contact with protein amino acids including W134 and W213 (Figure ). Despite these structural similarities, TCE binds site s1 with a higher affinity. Given TCE’s favorable interactions with lipids (Figure ), the predominant hydrophobic nature of site s1, as highlighted in Figure , makes sense of the result.

Figure 7

Characterization of sites s1 and s2. Per-site amino acid contacts (a) and number of water molecules (b) within 3 Å of bound TCE. Data was computed over independent equilibrium MD simulations (30 ns) of the ligand–protein bound state.

According to eq , we also calculated the standard binding free energy corresponding to both simultaneously bounded states s1 and s2 (unique saturation): ΔGs1,s2° = −5.45 kcal/mol. This value, comparable to the standard binding free energy for POPC (ΔGPOPC° = −4.32 ± 0.69 kcal/mol), evidences TCE’s greater affinity for the receptor rather than for the membrane. By these results, TCE is expected to bind protein targets with higher affinities, opening the possibility that it may bind ion channels when partitioning the membrane.

Conclusions

A fine atomistic CHARMM36-compatible model for TCE, a sedative hypnotic drug, is presented in this article. Model development targets gas-phase conformations and molecular electrostatic potential with individual TIP3P water molecules via weak hydrogen bonding. Validation is ensured, as developed parameters appropriately reproduce ligand’s physical and chemical properties, including liquid bulk properties, lipid partitioning, and interaction to protein targets.

Overall, our results support that the presented TCE force-field parameters are robust and likely to be useful in a large series of in silico biological studies. TCE is known experimentally to interact favorably with albumin and phospholipids as judged respectively from its dissociation constant and water–oil partition coefficient. Figures and 5 not only make sense of these experimental facts but also provide a detailed molecular model for ligand interaction to these substrates. Such a model is particularly useful to interpret and design novel experiments for characterization of ligand binding to membrane and albumin with potential consequences for our understanding of TCE action in membrane-embedded proteins. In this regard, our contribution opens the possibility to explore novel problems regarding the interaction of small ligands and ion channels, a field of immense interest.[22] Particular attention might be driven to K2P channels as the scientific literature still lacks information about their modulation mechanism by TCE and related molecules.

Computational Methods

Parametrization Adjustment

The MATCH atom-typing toolset was used to set initial atom parameters for TCE.[23] The optimization protocol proceeded in stages: acquisition of quantum mechanics (QM) and molecular mechanics (MM) data, comparison of conformational properties, and refinements. Distributed as a VMD plugin,[19] the force field toolkit[20] contains helpful scoring algorithms to fit MM to QM data[24] that allowed determination of accurate TCE parameters. In detail, the TCE equilibrium geometry at the QM level was defined as the molecular target. Next, atomic charge distributions were determined on the basis of QM interactions with individual TIP3P water molecules. For each donor or acceptor hydrogen bond, a typical linear geometry was built and then bond distances were optimized keeping fixed all other degrees of freedom. Ligand–water interaction energies were scaled by a factor of 1.16, whereas hydrogen bond lengths were shifted by an offset of −0.2 Å to yield appropriate bulk phase parameters.[25,26] QM partial atomic charges for the equilibrium geometry were used as a trial set for the charge optimization procedure. The objective function was then optimized on the basis of these QM interaction energies for TCE–water complexes until convergence was reached. During this iterative process, the QM dipole was allowed to be overestimated by 20–50% to reach consistency with the bulk phase.[27]

Bond and valence angle parameters were optimized by computing the energetic perturbation of small distortions from the equilibrium QM geometry along redundant internal coordinates. Optimization proceeded until convergence of bond and angle objective functions. Torsional parameters were obtained by fitting QM potential energy surfaces (PESs) for TCE dihedrals Cl–C2–C1–O and H–O–C1–C2. Torsion scans at the QM level were conducted by rotating dihedrals from −180 to +180° with a step size of 5°, providing an energy function in the form

where V is the fitted barrier height (force constant), and k and δ are dihedral multiplicity and phase, respectively. Parameters were adjusted to fit QM PES until convergence of the root-mean-square deviation. Finally, Lennard-Jones (LJ) parameters used to treat van der Waals interactions were assigned by transferability of similar chemical types available at the CHARMM general force field (CGenFF). After running a full optimization cycle, the geometry was minimized via a 1000-step gradient to compare the developed model to the QM target. As a standard practice, TCE–water interactions were reoptimized using this MM geometry, whereas Hessian and PES remained unmodified and a new round of parametrization was performed to ensure self-consistency. Further rounds were not required because no significant improvement was obtained. All MM calculations were performed using NAMD 2.10.[28]

QM Calculations

Molecular geometry was optimized at the MP2/6-31G(d) level, and initial partial atomic charges were derived from MP2/6-31G(d) Merz–Kollman charges.[29] Water interaction profiles were optimized at the HF/6-31G(d) level. Although higher levels of QM theory may lead to more accurate geometries and hydrogen bond energies, the chosen level of theory for water interaction profiles maintains consistency with the CHARMM additive force field.[14,27] The QM Hessian matrix and the relaxed PES were also calculated at the MP2/6-31G(d) level. Constraints were imposed to the scanned dihedrals during molecule minimization. All electronic structure calculations were performed in Gaussian 09.[15]

MD Simulations

Condensed phase simulations were performed in an NPT ensemble at 298.15 K and 1 bar using Langevin dynamics and the Langevin piston algorithm as implemented in NAMD 2.10.[28] Periodic boundary conditions were applied to all MD simulations. The particle mesh Ewald method was employed to evaluate full electrostatics using a real-space grid spacing of 1.2 Å or less. The multiple-time-step r-RESPA integrator was used, with a base time step of 2 fs for short-range nonbonded forces and extended time of 4 fs for soft, long-range interactions. Pairwise nonbonded interactions were truncated at a distance of 12 Å with a smooth switching function above (10 Å).

The starting configuration for liquid phase simulations consisted of 216 TCE molecules placed with random orientation at the grid points of a cubic lattice (6.5 × 6.5 × 6.5 Å3). A 10 000-step minimization followed by gradual heating during 10 ps was applied to the whole system for equilibration. Equilibrium data was collected for 0.4 ns to compute the average periodic cell volume to estimate TCE density.

Under the ideal-gas assumption and negligible mechanical work in the liquid phase,[14] the enthalpy of vaporization at a given temperature T can be written aswhere and ⟨Ugas⟩ correspond, respectively, to the average potential energy of the molecule in liquid and gas phases. First, was estimated on the basis of the liquid phase simulation of N = 216 molecules of TCE. Then, each molecule’s final configuration in the liquid phase was isolated and subsequently simulated in gas phase. These individual simulations comprised 50 ps of equilibration with a friction coefficient of 5 ps–1 followed by 50 ps of data collection. Finally, ⟨Ugas⟩ was obtained by averaging over the potential energy relative to N = 216 configurations.

TCE solvation free energies in bulk water and oil phase were evaluated using TIP3P water and hexane molecules, respectively.[30] Both solvent models are available in CHARMM36. Solvation free energies were computed via free energy perturbation calculations (FEP).[31] In these alchemical transformations, the solute is decoupled from the environment by turning its intermolecular interactions off, using a scalar parameter λ. A soft-core potential was adopted to scale the nonbonded pair potential of the perturbed system according to a dual coupling parameter (van der Waals and electrostatics).[32] Alchemical transformations were split into 100 “windows” and carried out explicitly in both directions, i.e., decoupling (integration over λ from 0 to 1) and recoupling simulations (integration over λ from 1 to 0). Each window contains 4 ps of relaxation followed by 66 ps of data collection. The soft-core potential shift distance was set to 7.0 Å2 for water and 5.5 Å2 for hexane, after several optimizations by trial and error. In total, five independent simulations (replicas) were performed with different initial velocities and seeds for the stochastic term in the Langevin thermostat. Statistical errors related to solvation free energies were computed with the simple overlap sampling (SOS) algorithm.[18]

Partitioning to Lipid Bilayer

A TCE molecule was inserted in a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer model, available in CHARMM36. A lamella with dimensions of approximately 45.0 × 43.0 × 89.0 Å3 was built using 59 lipids. The POPC bilayer was hydrated by 3171 TIP3P water molecules. Free energy profile or potential of mean force (PMF) for TCE to cross the membrane was computed using the adaptive biasing force algorithm (ABF) as implemented in NAMD.[33,34] In the method, the reaction coordinate to be sampled ξ is discretized in i bins. PMF is then derived from a biasing force added to the system’s equations of motion to counteract average forces acting on each bin. Here, the reaction coordinate ξ, with bin size 0.1 Å, and defined as the distance z between TCE and the center of POPC lamella, was sampled over 300 ns simulation. The accumulated statistical error throughout the PMF was calculated according to the equation[35]where i and i are bin indices delimiting the ξ-interval [a, b]; Δt is the simulation time step; and, for each bin, n, τ, and ⟨ΔFξ2⟩, respectively, stand for the number of samples accrued, autocorrelation time, and variance of the instantaneous force Fξ.

The partition coefficient K was computed from the PMF following equation[36]where |a – b| is the system width from bulk to the membrane center, along z. Partition coefficient statistical uncertainty was estimated by propagating the average bin error Err[ΔG]/i.

Binding Affinity to BSA

Fluorescence quenching experiments suggest that TCE binds to bovine serum albumin (BSA) at two distinct sites, namely, in the vicinity of tryptophans 134 and 213.[13] To further characterize and validate the TCE model, molecular docking and free energy perturbation (FEP) calculations of TCE against BSA were performed and compared with experimental data.

AutoDock Vina software was used to dock the ligand into an ensemble of BSA equilibrium structures to allow for sampling of the protein’s configurational space.[37] Binding sites were delimited from the ensemble of docking poses in such a way that each binding site was defined as the effective volume outlined by protein residues with the largest number of contacts to docking poses. In this manner, site s1 is defined by residues L197, W213, S343, L346, S453, L480, and V481 and site s2 by residues K20, V40, K131, W134, and G135 (Figure ). Although docking was performed in vacuum, subsequent equilibrium and FEP simulations were conducted in the presence of explicit all-atom water molecules, thus taking environmental effects into account. Specifically, the simulated system comprised the protein receptor (BSA) embedded in a homogeneous reservoir, with diluted ligand (TCE) concentration. The protein is assumed to be in a well-defined conformational state in which it provides two distinct TCE binding sites, s1 and s2. Hence, the equilibrium constant for the process of bringing the ligand from the bulk into the bound state (s1, s2) can be solved by means of an MD/FEP approach.[38] The method requires a harmonic potential coupled to the ligandto guarantee that it is restrained to occupy the effective volume centered at the receptor binding site R*; the restraint k exerts no force when the ligand center of mass is within a distance R from the equilibrium position and exerts a harmonic restoring force when the ligand center of mass is out of this range. In our case, force constants were estimated from the root-mean-square fluctuation of docking solutions. The equilibrium binding constant for each individual site j, j ∈ (s1, s2), can be written asBy definition, W*(R) is the reversible work to transfer a single ligand from the gas phase to the respective bound state and μ̅ refers to its solvation free energy. Thereon, assuming that the binding sites are independent, an aggregate affinity constant, for both binding sites to be simultaneously occupied by a single ligand each, can be defined as the product of the individual constants[38]From eq , the standard binding free energy relative to bringing a single ligand from a reference standard reservoir concentration into the respective target site j can be defined aswhere C° = 1 M or in units of number density C° = (1660 Å3)−1. Finally, from eq , we can derive the standard binding free energy relative to bringing two indistinguishable ligands from a reference standard reservoir concentration into simultaneously occupied binding sites s1 and s2To compute binding affinities described above, alchemical transformations were conducted by decoupling and recoupling the TCE gauche conformer from each BSA binding site separately. Starting from protein–TCE equilibrated systems as resolved from docking, decoupling/recoupling was carried out by varying the λ coupling parameter in steps of 0.01, amounting to 100 windows per transformation. Transformations in each direction spanned 2 ps of relaxation followed by 6.4 ns of data collection, thus totaling 12.8 ns of collection per replica, per site. For the purpose of improving statistics, site-specific FEP estimates and the associated statistical errors were determined from five independent decoupling/recoupling runs. Restraints on the ligand’s center of the mass were included in the calculations to ensure appropriate sampling within the binding site so that TCE’s chemical potential remains well defined during the final decoupling or initial recoupling stages.[39,40]