Hbond Autocorrelation functions and Lifetime of protein backbone water interaction

[Asakerei]

Hello!

I am currently investigating the effects co-solvent on H-bond donation/acceptance properties of water as it pertains to protein stability. I am interested in finding the H-bond lifetime (from ACF) between backbone oxygen and water (solvent donation) and backbone nitrogen and water (acceptance). I am having difficulty obtaining reproducible time constant at different trajectory steps of my simulation (each step is at least 3ns apart [1-9.trr]). My script can be found bellow (solvent donation).

any suggestions as to what I can change in the script?

Thanks in advance!
IA

import MDAnalysis as mda
from MDAnalysis.analysis.hydrogenbonds import hbond_autocorrel
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Define input files and output naming pattern
tpr_file = "md-300.tpr"
trr_files = [f"{i}.trr" for i in range(1, 10)]  # 1.trr to 9.trr
output_prefix = "HN-water_"

for idx, trr in enumerate(trr_files, start=1):
    print(f"🔄 Processing {trr}...")

    # Load trajectory
    u = mda.Universe(tpr_file, trr)
    
    # Select atoms
    H = u.select_atoms('name HN N')
    O = u.select_atoms('name OW')
    N = u.select_atoms('name N or name HN')

    # Run Hydrogen Bond AutoCorrelation
    hb_ac = hbond_autocorrel.HydrogenBondAutoCorrel(
        u,
        acceptors=O,
        hydrogens=H,
        donors=H,
        bond_type='continuous',
        sample_time=50.0,  # Ensure long enough range
        nsamples=200,  # Reduce truncation
        nruns=100
    )

    print(f"✅ Running Hydrogen Bond AutoCorrelation on {trr}...")
    
    hb_ac.run()
    hb_ac.solve()

    # Extract results
    time = hb_ac.solution['time']
    results = hb_ac.solution['results']
    estimate = hb_ac.solution['estimate']

    # Save results to CSV
    csv_filename = f"{output_prefix}{idx}.csv"
    df = pd.DataFrame({"Time (ps)": time, "Results": results, "Estimate": estimate})
    df.to_csv(csv_filename, index=False)

    print(f"✅ Data saved: {csv_filename}")

    # Plot results
    plt.figure(figsize=(8, 5))
    plt.plot(time, results, 'ro', label="Results")
    plt.plot(time, estimate, label="Exponential Fit")
    plt.xlabel("Time (ps)")
    plt.ylabel("Hydrogen Bond AutoCorrelation")
    plt.legend()
    plt.title(f"Hydrogen Bond AutoCorrelation - {trr}")
    
    # Save plot
    plot_filename = f"{output_prefix}{idx}.png"
    plt.savefig(plot_filename, dpi=300)
    plt.close()
    
    print(f"📊 Plot saved: {plot_filename}")

print("✅ All trajectories processed successfully!")

[orbeckst]

I don't have an answer to your question.

Just the advice to start with a system where you know the answer, where you can reproduce what you expect/know to be true, and then apply your tools to unknown (and therefore interesting) systems.

Any answer may be highly system dependent and given that you are the expert on the system that you're simulating, you should be in the best position to answer. As a generic answer: Vary your various parameters, see if it's a sampling issue by using twice as much (or ten times as much) data.

(But nothing beats starting from a well-sampled system where you can get reproducible results and where you can then vary the parameters to better understand how sensitive your results are to sample_time, nsamples, ... and your trajectory length.)