Protein in water/glycerol

The following examples consider a system composed a protein solvated by a mixture of water and glycerol, built with Packmol. The simulations were performed with NAMD with periodic boundary conditions and a NPT ensemble at room temperature and pressure. Molecular pictures were produced with VMD and plots were produced with Julia's Plots library.

Image of the system of the example: a protein solvated by a mixture of glycerol (green) and water, at a concentration of 50%vv.

Index

• Data, packages, and execution
• MDDF, KB integrals, and group contributions
• Preferential interactions and m-values
• 2D density map
• 3D density map

Data, packages, and execution

The files required to run this example are available at [this link], and are:

• system.pdb: The PDB file of the complete system.
• glyc50_traj.dcd: Trajectory file. This is a 1GB file, necessary for running from scratch the calculations.

To run the scripts, we suggest the following procedure:

1. Create a directory, for example example1.
2. Unzip and copy the required data files above to this directory.
3. Launch julia in that directory, activate the directory environment, and install the required packages. This is done by launching Julia and executing:

   import Pkg 
   Pkg.activate(".")
   Pkg.add(["ComplexMixtures", "PDBTools", "Plots", "LaTeXStrings", "EasyFit"])
   exit()

4. Copy the code of each script in to a file, and execute with:

   julia -t auto script.jl

   Alternatively (and perhaps preferably), copy line by line the content of the script into the Julia REPL, to follow each step of the calculation. For a more advanced Julia usage, we suggest the VSCode IDE with the Julia Language Support extension.

MDDF, KB integrals, and group contributions

Here we compute the minimum-distance distribution function, the Kirkwood-Buff integral, and the atomic contributions of the solvent to the density. This example illustrates the regular usage of ComplexMixtures, to compute the minimum distance distribution function, KB-integrals and group contributions.

# Activate environment in current directory
import Pkg;
Pkg.activate(".");

# Run this once, to install necessary packages:
# Pkg.add(["ComplexMixtures", "PDBTools", "Plots", "LaTeXStrings"])

# Load packages
using ComplexMixtures
using PDBTools
using Plots, Plots.Measures
using LaTeXStrings

# The complete trajectory file can be downloaded from (1Gb):
# https://www.dropbox.com/scl/fi/zfq4o21dkttobg2pqd41m/glyc50_traj.dcd?rlkey=el3k6t0fx6w5yiqktyx96gzg6&dl=0

# The example output file is available at:
# 
# Load PDB file of the system
atoms = read_pdb("./system.pdb")

# Select the protein and the GLYC molecules
protein = select(atoms, "protein")
glyc = select(atoms, "resname GLYC")

# Setup solute and solvent structures
solute = AtomSelection(protein, nmols=1)
solvent = AtomSelection(glyc, natomspermol=14)

# Path to the trajectory file
trajectory_file = "./glyc50_traj.dcd"

# Run mddf calculation, and save results
results = mddf(trajectory_file, solute, solvent, Options(cutoff=15.0))
save(results, "glyc50_results.json")
println("Results saved to glyc50_results.json")

#
# Produce plots
#
# Default options for plots 
Plots.default(
    fontfamily="Computer Modern",
    linewidth=2,
    framestyle=:box,
    label=nothing,
    grid=false
)

#
# The complete MDDF and the Kirkwood-Buff Integral
#
plot(layout=(1, 2))
# plot mddf
plot!(results.d, results.mddf,
    xlabel=L"r/\AA",
    ylabel="mddf",
    subplot=1
)
hline!([1], linestyle=:dash, linecolor=:gray, subplot=1)
# plot KB integral
plot!(results.d, [ 1e-3 * kbi(results) 1e-3 * kbi(results; correction=:none) ], #to L/mol
    xlabel=L"r/\AA",
    ylabel=L"G_{us}/\mathrm{L~mol^{-1}}",
    label=[ "KBI (with finite-size corrections)" "KBI (raw)" ],
    subplot=2
)
# size and margin
plot!(size=(800, 300), margin=4mm)
savefig("./mddf.png")
println("Created plot mddf.png")

#
# Atomic contributions to the MDDF
#
hydroxyls = ["O1", "O2", "O3", "H1", "H2", "H3"]
aliphatic = ["C1", "C2", "HA", "HB", "HC", "HD"]
hydr_contrib = contributions(results, SolventGroup(hydroxyls))
aliph_contrib = contributions(results, SolventGroup(aliphatic))

plot(results.d, results.mddf,
    xlabel=L"r/\AA",
    ylabel="mddf",
    size=(600, 400)
)
plot!(results.d, hydr_contrib, label="Hydroxyls")
plot!(results.d, aliph_contrib, label="Aliphatic chain")
hline!([1], linestyle=:dash, linecolor=:gray)
savefig("./mddf_atom_contrib.png")
println("Created plot mddf_atom_contrib.png")

Output

The code above will produce the following plots, which contain the minimum-distance distribution of glycerol relative to the protein, and the corresponding KB integral:

The plot shows the KBI computed with finite-size corrections (kbi(results)), which is the default output and the best estimate of the integral, and the uncorrected (raw - truncated) KBI. The uncorrected KBI can be useful to interpret how local interactions contribute to the exclusion or accumulation of the solvent at short distances. See the Kirkwood-Buff integrals section for further information on the corrections.

The distribution function can be decomposed into the contributions of the hydroxyl and aliphatic groups of glycerol:

options = Options(cutoff=15.0, stride=5)
mddf(trajectory_file, solute, solvent, options)

Preferential interactions and m-values

The KB integral of a single component of the solvent is difficult to interpret by itself: for a large solute, it is dominated by the excluded volume of the solute, and is thus large and negative for all components (about -20 L mol⁻¹ for both water and glycerol here). The thermodynamically relevant quantity is the difference between the KB integrals of the cosolvent and of water, which defines the preferential interaction parameter of the cosolvent ($c$, glycerol) relative to water ($w$), around the protein ($p$) (Shimizu, 2004; Martínez and Shimizu, 2017):

\[\Gamma_{pc} = \rho_c \left( G_{pc} - G_{pw} \right)\]

where $\rho_c$ is the bulk molar concentration of the cosolvent. $\Gamma_{pc} > 0$ indicates that the cosolvent accumulates on the protein surface relative to water (preferential solvation by the cosolvent), and $\Gamma_{pc} < 0$ that it is preferentially excluded. The excluded volume contributions cancel in the difference, and $\Gamma_{pc}$ can be interpreted as the excess number of cosolvent molecules in the solvation domain of the protein, relative to the composition of the bulk solution.

$\Gamma_{pc}$ is related to the dependence of the chemical potential of the protein, $\mu_p$, on the activity of the cosolvent, $a_c$, at constant temperature and pressure:

\[\left(\frac{\partial \mu_p}{\partial \ln a_c}\right)_{T,P} = -RT\,\Gamma_{pc}\]

and, therefore, on the concentration of the cosolvent (Shimizu, 2004):

\[\left(\frac{\partial \mu_p}{\partial \rho_c}\right)_{T,P} = -\frac{RT \left( G_{pc} - G_{pw} \right)}{1 + \rho_c\left( G_{cc} - G_{cw} \right)}\]

The denominator, $\left(\partial \ln a_c / \partial \ln \rho_c\right)^{-1}$, accounts for the non-ideality of the water-cosolvent mixture, and depends on the KB integrals between the components of the solvent, $G_{cc}$ and $G_{cw}$. If the mixture is ideal, $G_{cc} = G_{cw}$, and the denominator is one.

The m-value of a process involving the protein (for example, unfolding) is the derivative of its free energy with respect to the cosolvent concentration, and is thus obtained from the difference of the quantities above in the two states (for example, unfolded, U, and folded, F, states):

\[m = \frac{\partial \Delta G_{F \to U}}{\partial \rho_c} = -\frac{RT \left[ \left( G_{pc} - G_{pw} \right)^U - \left( G_{pc} - G_{pw} \right)^F \right]}{1 + \rho_c\left( G_{cc} - G_{cw} \right)}\]

The simulation of a single state, as here, provides the solvation contribution of that state to the m-value.

Computing the preferential interaction parameter

To compute $\Gamma_{pc}$, the KB integral of water relative to the protein is required, in addition to that of glycerol, computed above. The following script computes the MDDF of water from the same trajectory, and computes $\Gamma_{pc}$ and the derivative of the chemical potential of the protein, assuming that the water-glycerol mixture is ideal:

# Activate environment in current directory
import Pkg;
Pkg.activate(".");

# Run this once, to install necessary packages:
# Pkg.add(["ComplexMixtures", "PDBTools", "Plots", "LaTeXStrings"])

# Load packages
using ComplexMixtures
using PDBTools
using Plots, Plots.Measures
using LaTeXStrings
using Statistics

# Load PDB file of the system
atoms = read_pdb("./system.pdb")

# The results for glycerol were computed in script1.jl
glyc_results = load("./glyc50_results.json")

# Compute the MDDF of water relative to the protein, from the same trajectory
protein = select(atoms, "protein")
water = select(atoms, "water")
solute = AtomSelection(protein, nmols=1)
solvent = AtomSelection(water, natomspermol=3)
trajectory_file = "./glyc50_traj.dcd"
water_results = mddf(trajectory_file, solute, solvent, Options(cutoff=15.0))
save(water_results, "water_results.json")
println("Results saved to water_results.json")

#
# KB integrals (converted from cm³/mol to L/mol), and the bulk concentration of glycerol (mol/L)
#
d = glyc_results.d
G_pc = glyc_results.kb / 1000 # protein-glycerol
G_pw = water_results.kb / 1000 # protein-water
ρ_c = overview(glyc_results).density.solvent_bulk

# Preferential interaction parameter, as a function of the distance
Γ = ρ_c .* (G_pc .- G_pw)

# Converged values: average in the region where the KB integrals are stable
converged = findall(r -> 10.0 <= r <= 15.0, d)
Γ_conv = mean(Γ[converged])
ΔG_conv = mean(G_pc[converged] .- G_pw[converged])

# Derivative of the chemical potential of the protein with respect to the glycerol
# concentration, assuming that the water-glycerol mixture is ideal (kcal mol⁻¹ M⁻¹)
T = 298.15 # K
RT = 1.987204e-3 * T # kcal/mol
dμ_dρ_ideal = -RT * ΔG_conv

println("Bulk concentration of glycerol: $(round(ρ_c; digits=2)) mol/L")
println("G_pc - G_pw = $(round(ΔG_conv; digits=1)) L/mol")
println("Preferential interaction parameter, Γ = $(round(Γ_conv; digits=1))")
println("∂μ_p/∂ln(a_c) = -RTΓ = $(round(-RT * Γ_conv; digits=1)) kcal/mol")
println("∂μ_p/∂ρ_c (ideal mixture) = $(round(dμ_dρ_ideal; digits=1)) kcal mol⁻¹ M⁻¹")

#
# Plots
#
Plots.default(
    fontfamily="Computer Modern",
    linewidth=2,
    framestyle=:box,
    label=nothing,
    grid=false
)
plot(layout=(1, 2))
plot!(d, G_pc, label="Glycerol", subplot=1)
plot!(d, G_pw, label="Water", subplot=1)
plot!(xlabel=L"r/\AA", ylabel=L"G_{pi}/\mathrm{L~mol^{-1}}", subplot=1)
plot!(d, Γ, subplot=2)
hline!([0], linestyle=:dash, linecolor=:gray, subplot=2)
plot!(xlabel=L"r/\AA", ylabel=L"\Gamma_{pc}", subplot=2)
plot!(size=(800, 300), margin=4mm)
savefig("./preferential_interaction.png")
println("Created plot preferential_interaction.png")

Output

The KB integrals of glycerol and water relative to the protein, and the corresponding preferential interaction parameter, as a function of the distance to the protein, are:

The KB integrals, and thus $\Gamma_{pc}$, are stable for distances greater than about 8 Å, indicating that the bulk solution is reached at these distances from the protein surface. The script prints the values averaged between 10 and 15 Å:

Bulk concentration of glycerol: 5.95 mol/L
G_pc - G_pw = 2.1 L/mol
Preferential interaction parameter, Γ = 12.6
∂μ_p/∂ln(a_c) = -RTΓ = -7.4 kcal/mol
∂μ_p/∂ρ_c (ideal mixture) = -1.3 kcal mol⁻¹ M⁻¹

In this simulation, thus, glycerol is preferentially bound to the protein: there are about 13 more glycerol molecules in the solvation domain of the protein than expected from the bulk composition. The chemical potential of the protein decreases with the increase of the glycerol concentration (assuming the ideality of the water-glycerol mixture). The negative values of $\Gamma_{pc}$ at short distances reflect the fact that water, being smaller, approaches the protein surface more closely than glycerol. The 2D and 3D density maps below illustrate which regions of the protein are responsible for the interactions with glycerol.

Accounting for the non-ideality of the solvent

The non-ideality factor, $1 + \rho_c(G_{cc} - G_{cw})$, can be obtained from experimental activity data of the water-cosolvent mixture, or from a simulation of the water-cosolvent mixture without the protein, at the same composition. In this case, the solvent components are small molecules, and each can be represented by a single atom (for example, the central carbon of glycerol, C2, and the oxygen atom of water, OH2). Thus, the KB integrals are computed from radial distribution functions, and are obtained with the kbi function. For instance, assuming that mixture.pdb and mixture.dcd are the files of a simulation of the water-glycerol mixture:

using ComplexMixtures, PDBTools
atoms = read_pdb("mixture.pdb")
glyc_C2 = AtomSelection(select(atoms, "resname GLYC and name C2"), natomspermol=1)
water_O = AtomSelection(select(atoms, "water and name OH2"), natomspermol=1)
# glycerol-glycerol: solute and solvent are the same
R_cc = mddf("mixture.dcd", glyc_C2, Options(cutoff=20.0))
# glycerol-water 
R_cw = mddf("mixture.dcd", glyc_C2, water_O, Options(cutoff=20.0))
# KB integrals (L/mol), and the bulk concentration of glycerol (mol/L)
G_cc = kbi(R_cc) / 1000
G_cw = kbi(R_cw) / 1000
ρ_c = overview(R_cc).density.solvent_bulk
# Non-ideality factor, as a function of L (take the value where it is stable)
nonideality = @. 1 + ρ_c * (G_cc - G_cw)

The value of $\partial \mu_p / \partial \rho_c$ computed above, assuming ideality, must then be divided by this factor.

2D density map

In this followup from the example above, we compute group contributions of the solute (the protein) to the MDDFs, split into the contributions each protein residue. This allows the observation of the penetration of the solvent on the structure, and the strength of the interaction of the solvent, or cosolvent, with each type of residue in the structure. The ResidueContributions and Plots.contourf auxiliary functions, documented here, are used:

# Activate environment in current directory
import Pkg;
Pkg.activate(".");

# Run this once, to install necessary packages:
# Pkg.add(["ComplexMixtures", "PDBTools", "Plots", "LaTeXStrings"])

# Load packages
using ComplexMixtures
using PDBTools
using Plots

# Load PDB file of the system
atoms = read_pdb("./system.pdb")

# Select the protein and the GLYC molecules
protein = select(atoms, "protein")

# Load example output file (computed in the previous script)
example_output = "./glyc50_results.json"
results = load(example_output)

#
# Plot a 2D map showing the contributions of some residues
#
residue_contributions = ResidueContributions(
  results, 
  select(protein, "resnum >= 70 and resnum <= 110")
)
contourf(residue_contributions; oneletter=true)
savefig("./density2D.png")
println("Created plot density2D.png")

Output

The code above will produce the following plot, which contains, for each residue, the contributions of each residue to the distribution function of glycerol, within 1.5 to 3.5 $\mathrm{\AA}$ of the surface of the protein.

3D density map

In this example we compute three-dimensional representations of the density map of Glycerol in the vicinity of a set of residues of a protein, from the minimum-distance distribution function. For further information about the grid3D function, go to the 3D grid map section of the manual.

import Pkg;
Pkg.activate(".");
using PDBTools
using ComplexMixtures

# PDB file of the system simulated
atoms = read_pdb("./system.pdb")
protein = select(atoms, "protein")

# Load results of a ComplexMixtures run
results = load("./glyc50_results.json")

# Compute the 3D density grid. Here we use dmax=3.5 such that the
# grid is not too large. Provide a file name (e.g. "./grid.pdb") as
# the third argument to also write the grid to a PDB file.
grid = grid3D(results, atoms; dmin=1.5, dmax=3.5)

# Regions associated with the hydrogen-bonding peak of the MDDF (d < 2.0 Å),
# with at least 10% of the maximum contribution. The protein is shown as a 
# white cartoon, the grid points as dots colored by their relative contribution
# from white (zero) to red (≥ 0.4), and the GLU and ASP residues as green balls and sticks
hbonds = filter(p -> occup(p) < 2.0 && beta(p) > 0.1, grid)
view_hbonds = visualize(
    protein => (color="white",),
    hbonds => (style=:dots, color_by=beta.(hbonds), colormap=:rwb, color_range=(0.4, -0.4)),
    protein => (selection="resname GLU ASP", style=:ballandstick, color="green"),
)
save("./grid_hbonds.html", view_hbonds)

# Regions associated with the second peak of the MDDF (2.0 Å < d < 3.5 Å),
# with the THR, ASN, and TRP residues as balls and sticks
second_peak = filter(p -> 2.0 < occup(p) < 3.5 && beta(p) > 0.1, grid)
view_second_peak = visualize(
    protein => (color="white",),
    second_peak => (style=:dots, color_by=beta.(second_peak), colormap=:rwb, color_range=(0.4, -0.4)),
    protein => (selection="resname THR ASN TRP", style=:ballandstick),
)
save("./grid_second_peak.html", view_second_peak)

Here, the MDDF is decomposed at each distance according to the contributions of each solute (the protein) atom. The grid is created such that, at each point in space around the protein, it is possible to identify:

1. Which atom is the closest atom of the solute to that point.

2. Which is the contribution of that atom (or residue) to the distribution function.

Therefore, by filtering the 3D density map at each distance one can visualize over the solute structure which are the regions that mostly interact with the solvent of choice at each distance. The views are created with the visualize function of PDBTools (version 3.41.0 or greater), which renders them with 3Dmol.js (Rego and Koes, 2015), and the script saves them as the grid_hbonds.html and grid_second_peak.html files, which can be opened in any web browser. Typical views are shown below (rotate, zoom, and hover over the points to identify the protein atoms):

In the view above, the points in space around the protein are selected with the following properties: distance from the protein smaller than 2.0Å and relative contribution to the MDDF at the corresponding distance of at least 10% of the maximum contribution. Thus, we are selecting the regions of the protein corresponding to the most stable hydrogen-bonding interactions. The protein is shown as a white cartoon, the Glutamic and Aspartic acid residues as green balls and sticks, and the grid points as dots, colored by the contribution to the MDDF, from white to red. Thus, the most reddish points correspond to the regions where the most stable hydrogen bonds were formed.

Hovering over those points we obtain which are the atoms of the protein contributing to the MDDF at that region (for example, ASP14:P OD1 for the OD1 atom of residue 14). In particular, the strongest red regions correspond to the carboxylate groups of Aspartic and Glutamic acid residues.

The view below displays the most important contributions to the second peak of the distribution, corresponding to distances from the protein between 2.0 and 3.5Å, with the Threonine, Asparagine, and Tryptophan residues shown as balls and sticks:

Notably, the regions involved are different from the ones forming hydrogen bonds, indicating that non-specific interactions with the protein (and not a second solvation shell) are responsible for the second peak.

How to run this example:

Assuming that the input files are available in the script directory, just run the script with:

julia -t auto script3.jl

Alternatively, open Julia and copy/paste the commands of script3.jl or use include("./script3.jl"). These options will allow you to remain on the Julia session with access to the grid data structure that was generated, and to display the views interactively (in VSCode or in the browser, if running from the REPL).

Visualization with VMD

The grid can also be written to a PDB file and visualized with other software. An example of the visualization of the grid of this example with VMD is available here.