Site occupancy

The occupancy function computes, for each frame of a simulation, which solvent molecules are found within a given cutoff distance of a binding site, defined by a set of atoms (for example, the atoms lining a cavity, or the surface of a protein).

The distance considered is the minimum distance between the site atoms and the atoms of each solvent molecule, computed with the same minimum_distances! machinery used by coordination_number.

The function returns an Occupancy object, wrapping, for each frame, the list of solvent molecules found at the site, and the corresponding site-solvent distances. This object can then be used to compute simple summary statistics (with mean), to characterize how long, on average, individual solvent molecules remain at the site (with intermittent_correlation), or to resolve that residence time as a function of the distance to the site (with intermittent_correlation_profile and residence_time).

Example: TMAO occupancy at a protein binding site

Here we use the MolSimToolkit.Testing test data, a short NAMD trajectory of a protein solvated by water and TMAO, and compute the occupancy of the protein surface by TMAO molecules, using a 3 Å cutoff:

using MolSimToolkit, PDBTools, MolSimToolkit.Testing, Plots

sim = Simulation(Testing.namd2_pdb, Testing.namd2_traj)
protein = select(get_atoms(sim), "protein")
tmao = select(get_atoms(sim), "resname TMAO")

occ = occupancy(
    sim, protein, tmao;
    solvent_natomspermol=14, cutoff=3.0,
    show_progress=false,
)

-------------------------------------------------------------------
Occupancy data:
-------------------------------------------------------------------
Total number of solvent molecules: 181
Maximum distance to the site (dmax): 3.0
Mean occupancy: 5.3
Minimum occupancy: 2
Maximum occupancy: 10
-------------------------------------------------------------------

The list field of the Occupancy object contains, for each frame, the indices of the TMAO molecules (1-based, relative to the tmao selection) found within the cutoff distance of the protein:

occ.list

20-element Vector{Vector{Int64}}:
 [16, 59, 63, 79, 81, 144, 177]
 [16, 175, 176]
 [16, 53, 106, 131]
 ⋮
 [104, 179]
 [70, 77, 97, 125, 128, 168]

The distances field carries the corresponding minimum distances between the protein and each of those molecules, such that occ.distances[i][j] is the distance of the molecule occ.list[i][j]:

occ.distances

20-element Vector{Vector{Float64}}:
 [2.193289015874046, 2.1708631961351528, 2.1099792336190184, 2.6644466961410926, 2.243952690002156, 2.3923658051170125, 1.6032937405182721]
 [1.9506459897129325, 2.8333266769863195, 2.899149006265755]
 [1.8145360451686594, 2.810336352428233, 2.6114751082310885, 1.9895821730783263]
 ⋮
 [2.0022412140756876, 2.19089352500517]
 [2.9576917947458146, 1.769920350132631, 2.1694941674417874, 2.2090029171689043, 2.1467604947077152, 2.343788354119454]

n_solvent_molecules is the total number of TMAO molecules considered, and dmax is the cutoff used, that is, the maximum distance at which a molecule is considered to be at the site:

occ.n_solvent_molecules, occ.dmax

(181, 3.0)

Average occupancy

The average number of TMAO molecules found at the protein surface, per frame, is obtained with mean:

mean(occ)

5.3

Intermittent correlation of the occupancy

The intermittent_correlation function, applied to an Occupancy object, estimates the probability that a TMAO molecule found at the site at some frame is still found there delta frames later. This probability decays as molecules exchange between the site and the bulk solvent, and can be used to characterize the residence time of the solvent at the site:

ic = intermittent_correlation(occ; maxdelta=4, show_progress=false)

5-element OffsetArray(::Vector{Float64}, 0:4) with eltype Float64 with indices 0:4:
 1.0
 0.39
 0.23469387755102042
 0.12903225806451613
 0.033707865168539325

plot(MolSimStyle,
    0:4, parent(ic), # parent(c) is required for c is an OffsetArray.
    xlabel="Delta (frames)", ylabel="Probability",
    linewidth=2, marker=:circle,
    label="TMAO at protein surface",
)

Characteristic residence time

In this case, the decay of the intermittent correlation function can be fit to a double exponential, $c(\delta) = a\exp(-\delta/\tau)$, using EasyFit.fitexpdecay (already a dependency of MolSimToolkit.jl), to extract a characteristic residence time $\tau$, in units of frames:

using EasyFit
using JuMP, Ipopt # required for constrained non-linear fitting

fit = fitexpdecay(ic; n=2, c=mean(occ)/occ.n_solvent_molecules)
tau = fit.b

2-element Vector{Float64}:
 0.05556530042171695
 1.432360418011985

The constant term is set to the ratio of the mean coordination number and the total number of solvent molecules (correlation at long times).

plot(MolSimStyle,
    0:4, parent(ic),
    seriestype=:scatter,
    xlabel="Delta (frames)", ylabel="Probability",
    marker=:circle,
    label="TMAO at protein surface",
)
plot!(fit.x, fit.y, linewidth=2, label="Bi-Exponential fit")

Here the trajectory is very short (20 frames), so this characteristic time is shown only for illustration: with longer, production-quality trajectories, $\tau$ provides a quantitative estimate of how long a solvent molecule typically remains bound to the site.

Residence time as a function of the distance

The correlation function above mixes molecules that are tightly bound to the protein with molecules that are only loosely associated with it, at the edge of the cutoff. Since the Occupancy object also stores the site-solvent distances, the correlation function can be resolved by distance, with intermittent_correlation_profile.

The distances from 0 to dmax are split into bins of constant width delta_r, and the lower edges of consecutive bins are displaced by step_r. Using step_r smaller than delta_r produces overlapping bins, and thus a smooth (quasi-continuous) profile. For each bin, the correlation at delta is the probability that a molecule found in that bin at some frame is found in the same bin delta frames later.

Here we use a larger cutoff (6 Å), such that the profile covers the first solvation shells of the protein:

occ6 = occupancy(
    sim, protein, tmao;
    solvent_natomspermol=14, cutoff=6.0,
    show_progress=false,
)

profile = intermittent_correlation_profile(occ6;
    delta_r=1.0, step_r=0.25, maxdelta=4,
    show_progress=false,
)

-------------------------------------------------------------------
Intermittent correlation profile:
-------------------------------------------------------------------
Number of distance bins: 21
Bin width (delta_r): 1.0
Bin step (step_r): 0.25
Maximum distance (dmax): 6.0
Range of bin centers: 0.5 - 5.5
Maximum delta (frames): 4
Observations per bin: minimum = 0, maximum = 89
-------------------------------------------------------------------

The r field contains the center of each bin, correlations the corresponding correlation functions, and counts the number of observations found in each bin (bins with no observations produce NaN correlations):

plot(MolSimStyle,
    0:4, [parent(c) for c in profile.correlations[5:5:end]],
    xlabel="Delta (frames)", ylabel="Probability",
    linewidth=2,
    labels=hcat(["r = $(r) Å" for r in profile.r[5:5:end]]...),
)

The 50% residence time

The residence_time function returns, for each bin, the time at which the correlation function falls below a given threshold. With the default threshold=0.5, this is the time after which half of the molecules initially in the bin have left it: the 50% residence time. The value is obtained by linear interpolation between the two delta-steps that bracket the threshold, and is NaN for bins in which the correlation never falls below it (or for empty bins).

t50 = residence_time(profile; threshold=0.5)

plot(MolSimStyle,
    profile.r, t50,
    xlabel="Distance to the protein / Å",
    ylabel="50% residence time / frames",
    linewidth=2, marker=:circle,
    label=nothing,
)

Use the dt keyword to convert the result from frames to time units, for instance residence_time(profile; dt=0.1) if consecutive frames are separated by 0.1 ns.

As above, the trajectory used here is very short (20 frames), so the profile is shown only for illustration.

Reference functions

Occupancy

occupancy(
    sim::Simulation,
    site::AbstractVector{<:PDBTools.Atom},
    solvent::AbstractVector{<:PDBTools.Atom};
    solvent_natomspermol::Integer,
    cutoff::Real,
    show_progress::Bool = true,
)

julia> using MolSimToolkit, PDBTools, MolSimToolkit.Testing

julia> sim = Simulation(Testing.namd2_pdb, Testing.namd2_traj; frames=1:5);

julia> protein = select(get_atoms(sim), "protein");

julia> tmao = select(get_atoms(sim), "resname TMAO");

julia> occ = occupancy(sim, protein, tmao; solvent_natomspermol=14, cutoff=3.0, show_progress=false);

julia> length.(occ.list)
5-element Vector{Int64}:
 7
 3
 4
 5
 6

julia> occ.dmax
3.0

julia> occ.distances[2] # distances of the molecules of occ.list[2]
3-element Vector{Float64}:
 1.9506459897129325
 2.8333266769863195
 2.899149006265755

mean(occupancy::Occupancy)

julia> using MolSimToolkit, PDBTools, MolSimToolkit.Testing

julia> sim = Simulation(Testing.namd2_pdb, Testing.namd2_traj; frames=1:5);

julia> protein = select(get_atoms(sim), "protein");

julia> tmao = select(get_atoms(sim), "resname TMAO");

julia> occ = occupancy(sim, protein, tmao; solvent_natomspermol=14, cutoff=3.0, show_progress=false);

julia> mean(occ)
5.0