Seeing Stability at the Atomic Level

BsLipA (Bacillus subtilis Lipase A) is a 181-residue alpha/beta hydrolase used in detergent formulation, biodiesel production, and food processing. Its catalytic triad (Ser77, His156, Asp133) cleaves ester bonds. Like many industrial enzymes, its utility is limited by thermal stability.

The computational consensus flagged F58I as a top stabilizing candidate. The question this work answers: does F58I actually change the dynamic behavior of the protein under thermal stress? And does the active site geometry remain undistorted?

Static predictions score a single frozen structure. They cannot tell you how the protein breathes under heat, when conformational excursions begin, whether the active site geometry survives, or how reproducible the effect is across independent simulations.

Molecular dynamics answers all of these. It simulates every atom in the protein plus thousands of water molecules, computing forces and updating positions every 2 femtoseconds. Over 40 nanoseconds (20 million timesteps per run) we watch the protein vibrate, flex, and at high enough temperatures, begin to explore conformational space beyond the native basin.

Static stability prediction is like checking if a building's blueprints show enough steel beams. MD is like putting the building on a shake table and watching what happens.

Rather than simple geometric side-chain replacement, the platform used IntelliFold v2 and Chai-1 to cross-validate the Ile58 rotamer placement. Ten IntelliFold diffusion samples and 5 Chai-1 models (15 independent predictions) all converged on the same rotameric state, with highly consistent side-chain geometry (chi1 ≈ −70°, chi2 ≈ 170°).

The validated Ile58 side chain was then transplanted onto the wild-type crystal backbone (PDB 1ISP, 1.3 Å resolution). Both WT and F58I start from the identical crystal backbone, differing only at position 58. This eliminates starting-structure bias: any difference we observe is due to the mutation, not the starting structure quality.

AI structure prediction here served as rotamer cross-validation: it confirmed that the Ile58 placement was energetically consistent across two independent tools, giving confidence in the starting structure before committing GPU hours to simulation.

The simulation campaign used GROMACS 2024.4 with CUDA acceleration on NVIDIA A100 GPUs, the CHARMM36m force field with TIP3P water (the community standard for mutation-sensitive packing studies).

Total runs

28 independent simulations (1,120 ns aggregate)

300K (control)

5 WT + 8 F58I replicates

350K (thermal stress)

8 WT + 7 F58I replicates

Production MD per run

40 ns

Force field

CHARMM36m + TIP3P

Hardware

NVIDIA A100 GPUs

RMSD measures how far the protein backbone drifts from its starting structure over time. At 300K, both variants are stable with RMSD around 0.10 nm, with no significant difference. Under thermal stress at 350K, the gap appears: F58I shows 13% lower RMSD (0.130 vs 0.150 nm). The temperature scaling is notable: the RMSD gap grows 2.4× from 300K to 350K.

The variance across replicates reveals a suggestive pattern. Wild-type RMSD variance increases 8× from 300K to 350K. F58I variance barely changes.

At 350K, 6 of 8 wild-type replicates cluster tightly between 0.122-0.134 nm, while 2 show moderately elevated RMSD values of 0.207 and 0.225 nm, enhanced conformational fluctuations that are not observed in any of the 7 F58I replicates.

To be precise about what this means: RMSD values of ~0.2 nm represent conformational excursions beyond the native basin, not partial unfolding. Genuine loss of fold would produce RMSD values of 0.4 nm or higher. Near a melting transition, stochastic sampling of larger-amplitude motions in some replicates but not others is expected behavior, and the observation that only WT replicates show this is consistent with F58I having a higher energy barrier to conformational excursion. However, with 8 replicates per variant, the 2-vs-0 pattern is suggestive rather than statistically conclusive.

Per-residue RMSF (root mean square fluctuation) reveals where the protein becomes more or less flexible. F58I is 6.2% more rigid globally at 300K, with a pattern of selective rigidification: the hydrophobic core tightens while surface loops retain their natural flexibility.

Both variants show similar catalytic triad geometry (Ser77, His156, Asp133) at each temperature. F58I does not introduce obvious distortion of the active site. This is consistent with the structural predictions: all 15 independent models confirmed an intact triad in the static F58I structure.

Note: preserved geometry in simulation supports the absence of active-site disruption, but is not proof of catalytic activity; that requires experimental enzyme kinetics.

Five observations emerge from this validation:

1. F58I shows reduced conformational mobility under thermal stress: lower mean RMSD, tighter replicate clustering, and 6.2% lower per-residue flexibility at 300K. The effect scales with temperature.
2. The stabilization mechanism is improved core packing: isoleucine's branched aliphatic side chain packs more efficiently than phenylalanine's aromatic ring at position 58, rigidifying the surrounding structure. The precise energetic contribution merits further investigation.
3. The replicate variance pattern at 350K is suggestive: wild-type shows conformational excursions that F58I does not. Larger replicate ensembles and longer simulations would be needed to establish this with statistical confidence.
4. No obvious active-site distortion: the catalytic triad maintains similar geometry in both variants, consistent with the mutation being structurally benign at the active site.
5. AI structure prediction provided reliable rotamer cross-validation: 15 independent predictions converged on the same Ile58 placement, confirming the starting structure before committing to simulation.

This validation shows the pipeline and produces encouraging initial results. A publication-grade study would extend it in several directions:

• More replicates: 20+ per condition would provide the statistical power to quantify the conformational excursion frequency with confidence intervals.
• Longer trajectories: 100+ ns per run would improve convergence and reveal whether the elevated-RMSD replicates are transient fluctuations or progressive events.
• Per-residue packing analysis: quantifying the local contact density and van der Waals energy around position 58 to confirm the core packing mechanism.
• Predefined unfolding criteria: objective thresholds (Rg expansion, native contact loss, DSSP secondary structure) applied consistently across all replicates.
• Experimental validation: differential scanning fluorimetry or calorimetry to measure the actual ΔT_m, and enzyme kinetics to confirm catalytic activity is preserved.

The computational infrastructure, from stability prediction through structure validation to GPU-accelerated MD, is in place. Scaling it to publication depth is a matter of compute budget and collaboration scope, not capability.

This is what computational validation of a designed mutation looks like. A multi-tool consensus pipeline identifies a candidate. Structure prediction tools cross-validate the rotamer placement. Over a microsecond of GPU-accelerated molecular dynamics probes the dynamic consequences under realistic thermal conditions.

The results are encouraging: F58I shows reduced conformational mobility, higher intra-protein hydrogen-bond counts, and no obvious active-site distortion in this simulation set. The signal justifies deeper computational and experimental follow-up.

If you are engineering enzymes for thermal stability and want to validate your candidates computationally before committing to the lab, this is the kind of analysis the Crane Protein Intelligence Platform can run.

Let's talk→