LinB Sn2 mechnism

In this tutorial, we will model the first step of the catalytic reaction of a haloalkane dehalogenase, LinB, using a previously prepared molecular system. In this reaction, Asp108 acts as a nucleophile, attacking the carbon atom bonded to the chlorine atom of the substrate (SN2 reaction). This step, commonly referred to as the dechlorination step, leads to the formation of a covalent enzyme–substrate intermediate and the release of a chloride ion.

The workflow presented here combines three approaches:

1. Potential Energy Surface (PES) scanning to explore the reaction pathway.

2. Umbrella sampling to obtain a more statistically meaningful sampling of the reaction coordinate.

3. WHAM (Weighted Histogram Analysis Method) to reconstruct the corresponding free-energy profile.

This tutorial uses a deliberately small QC region to keep the computational cost low. The same workflow can be repeated using larger QC regions to investigate how the size of the quantum-mechanical region affects the calculated energy and free-energy profiles.

Step 1 – Opening the Prepared System

We will start from a molecular system that has already been prepared for the QM/MM calculation. To open the system, navigate to:

Main Menu → File → Open

Select the prepared YAML file:

LinB.yaml

After loading the system, inspect the molecular representation carefully.

The system already contains:

• A 15-atom QC region (AM1 level), displayed using a ball-and-stick representation.

• An outer layer of fixed atoms, represented in gray.

• The substrate carbon atoms highlighted in light purple, making them easier to distinguish from the chlorine atom.

• The chlorine atom represented using a darker green color.

Before proceeding, it is good practice to visually inspect the system and confirm that the QC region and the atoms involved in the reaction have been assigned correctly.

Step 2 – Modeling the Reaction Using a PES Scan

We will first explore the reaction pathway using a reaction-coordinate scan.

Open the reaction-coordinate scan tool through:

Simulation → Reaction Coordinate Scans

Defining the reaction coordinate:

Select the atoms involved in the reaction using the Picking mode.

The figure below illustrates the atom selection and the scan setup used in this calculation.

For this reaction, we will define the reaction coordinate as a linear combination of two distances:

RC = m₁ * d(OAsp108–C) − m₂ * d(C–Cl)

where m₁ and m₂ are two constants that depend on the masses of the selected atoms pk1 and pk3, respectively.  This type of reaction coordinate is particularly useful for describing a bond-breaking/bond-forming process because it simultaneously monitors the breaking C–Cl bond and the formation of the bond between Asp108 and the substrate carbon (Same approach to the one used in Tutorial 2, where we modeled a simple SN2 reaction).

Scan parameters

For this tutorial, 15 scan steps are sufficient to provide a quick overview of the reaction pathway.

Make sure to select a directory where you have write permission, since EasyHybrid will generate several output files during the calculation.

Once the setup is complete, click Run and wait for the scan to finish.

Tip: For a production-quality calculation, you may want to increase the number of scan points. Fifteen points are used here mainly to keep the tutorial fast and computationally inexpensive.

Step 3 –Importing and Inspecting the PES Scan Results

Once the scan has finished, we will import the resulting trajectory and inspect the structures generated along the reaction coordinate.

The first step is to import the trajectory into EasyHybrid.

Importing the trajectory

1. In the System Header, right-click to open the TreeView context menu.

2. Select Import Data...

3. Change the format to:

4. pkl folder – pDynamo trajectory

5. Select the folder generated by the reaction-coordinate scan.

6. If the calculation was completed successfully, EasyHybrid should automatically detect the log file associated with the trajectory.

7. Click Import.

The most important point at this stage is not the absolute energy, but whether the structures along the trajectory make chemical sense.

Check, for example, whether:

1. The Asp108 nucleophile approaches the substrate carbon.

2. The C–Cl bond becomes progressively longer.

3. The new Asp108–C bond becomes progressively shorter.

4. The chloride ion moves away from the substrate.

5. No unexpected structural distortions occur elsewhere in the system.

Step 4 – Analyzing the Potential Energy Profile

We can now inspect the potential energy profile generated by the scan.

Go to:

Analysis → PES Analysis

The PES Analysis window allows you to inspect the energy associated with each trajectory frame.

Use the scroll bar at the bottom of the window to move through the trajectory and inspect the corresponding molecular structures.

The resulting profile should provide a qualitative picture of the energetic changes associated with the dechlorination process.

For a relatively small system such as the one used in this tutorial, this is already a useful first look at the reaction pathway.

Important: A PES scan provides a potential-energy profile, not a free-energy profile. The next step will therefore use umbrella sampling to obtain a free-energy estimate along the same reaction coordinate.

Step 5 – Setting Up Umbrella Sampling

We will now calculate the free-energy profile using umbrella sampling. The structures generated by the PES scan will be used as the initial configurations for the umbrella-sampling windows. In this example, the scan contains 15 frames, so we will use 15 umbrella windows. Each frame provides the starting structure for one window. During the simulation, a harmonic biasing potential is applied to the reaction coordinate, keeping the system close to the target value associated with that window.

This allows the simulations to efficiently sample regions of the reaction coordinate that would otherwise be visited only rarely, particularly near high-energy regions.

Open the umbrella-sampling tool through:

Simulation → Umbrella Sampling

Alternatively, use the corresponding purple toolbar button.

Change the input type to:

• From trajectory (parallel)

This option allows the umbrella windows to be initialized directly from the trajectory generated by the PES scan.

Select the folder containing the trajectory frames generated during the scan.

Adjust the number of CPUs according to the computational resources available on your computer.

For this tutorial, we will use a simplified setup designed to keep the calculation short:

• Equilibration: 0 steps

• Production: 5000 steps

• Total sampling time: approximately 5 ps

The complete setup is shown in the figure below.

Note: No equilibration is performed here purely for tutorial purposes. For a production calculation, equilibration of each window is strongly recommended, and significantly longer sampling times should normally be used.

The umbrella-sampling calculation should take approximately 5–10 minutes on a modern multicore processor, although the actual time will depend on the QM method, number of CPUs, hardware, and system configuration.

Step 6 – Analyzing the Results with WHAM

Once all umbrella-sampling windows have finished, the free-energy profile can be reconstructed using the Weighted Histogram Analysis Method (WHAM).

Open the WHAM tool through:

Analysis → WHAM

Important WHAM Considerations

Before running WHAM, verify the following parameters.

1. Temperature

The temperature specified for the WHAM analysis must be identical to the temperature used in the molecular dynamics simulations.

Using inconsistent temperatures will result in an incorrect free-energy reconstruction.

2. Number of bins

The number of histogram bins should be sufficiently large to represent the sampled reaction coordinate accurately.

As a general guideline, the number of bins should be at least approximately twice the number of simulation windows.

In this tutorial, we have 15 windows, so 100 bins provides more than adequate resolution for this demonstration.

The optimal number of bins, however, depends on the amount of sampling and the distribution of the data.

3. Histogram overlap

A reliable WHAM free-energy profile requires sufficient overlap between neighboring histograms.

This is one of the most important aspects of umbrella sampling.

If neighboring windows are too far apart, or if the harmonic force constant is too large, each simulation may sample only a narrow region of the reaction coordinate. In that case, the histograms may have insufficient overlap, making the free-energy reconstruction unreliable.

Ideally, the histograms from neighboring windows should overlap substantially along the reaction coordinate.

Step 7 – Interpreting the Free-Energy Profile

After WHAM has completed, inspect the resulting free-energy profile.

The profile provides an estimate of the free-energy change along the selected reaction coordinate and can be used to identify:

• Reactant-like configurations.

• Transition-state-like regions.

• Product/intermediate-like configurations.

• The approximate free-energy barrier.

The relative free-energy difference between the initial and final states.

Remember that the free-energy profile obtained here is only as reliable as the sampling used to generate it. For this reason, the short simulations used in this tutorial should be considered demonstrative rather than production-quality.

Final Remarks

This tutorial illustrates a complete workflow for studying an enzymatic reaction using EasyHybrid:

Prepared QM/MM system

↓

Reaction Coordinate Scan

↓

Potential Energy Profile

↓

Umbrella Sampling

↓

WHAM

↓

Free-Energy Profile

It is important to emphasize that the QC region used in this tutorial was intentionally kept as small as possible to reduce the computational cost and make the tutorial accessible.

For a realistic computational study, it is worth repeating the calculation using larger QC regions and investigating whether the calculated reaction pathway, activation energy, and free-energy profile converge with respect to the size of the QM region.

You can therefore use this tutorial as a starting point and explore the effect of increasing the QC region. The computational cost will increase, but the resulting model may provide a more realistic description of the chemical environment involved in the catalytic reaction.

References:

Junjie Wang, Xiaowen Tang, Yanwei Li, Ruiming Zhang, Ledong Zhu, Jinfeng Chen, Yanhui Sun, Qingzhu Zhang, Wenxing Wang; Computational evidence for the degradation mechanism of haloalkane dehalogenase LinB and mutants of Leu248 to 1-chlorobutane. Phys. Chem. Chem. Phys. 2018; 20 (31): 20540–20547.