CosolvKit

The python package for creating cosolvent system.

The original paper is freely accessible at the link (https://pubs.acs.org/doi/10.1021/acs.jcim.4c01398).

Documentation

The installation instructions, documentation and tutorials can be found on http://cosolvkit.readthedocs.io/.

Installation

Cosolvkit package is available via conda and can be installed:

$ conda install --channel conda-forge::cosolvkit

or via mamba:

$ mamba install -c conda-forge cosolvkit

*Please note that Apple M1 chips are not supported by some of CosolvKit's dependencies. we recommend macOS users of Apple Silicon install the x86_64 version of MiniForge and run CosolvKit through Rosetta.

Finally, we can install the CosolvKit package via pip:

$ conda create -n cosolvkit -c conda-forge -f environment.yml
$ conda activate cosolvkit
$ pip instal cosolvkit

or directly download and install the source code from git: I highly recommend you to install the Anaconda distribution (https://www.continuum.io/downloads) if you want a clean python environnment with nearly all the prerequisites already installed. To install everything properly, you just have to do this (for faster installation, use mamba or micromamba instead of conda):

$ git clone https://github.com/forlilab/cosolvkit
$ cd cosolvkit
$ conda create -n cosolvkit -c conda-forge -f environment.yml
$ conda activate cosolvkit
$ pip install -e .

Example

To run an example of CosolvKit system creation run the command create_cosolvent_system -c cosolvkit/data/config.yaml. By default this example creates the cosolvent system without running any MD simulation or post processing analysis. The results can be found in cosolvkit/data/results.

Quick tutorial

The script create_cosolvent_system.py provides all the necessary tools to build a cosolvent system. The main entry point is the file config.yaml where all the necessary options are specified.

Argument	Type	Description	Default value	OPENMM	AMBER	GROMACS	CHARMM
cosolvents	list	List of cosolvent molecules. Each entry requires name, smiles, resname, and either concentration (mol/L) or copies (integer).	no default
md_engine	dict	Dict mapping engine name to list of forcefield XML files. Supported engines: [openmm, amber, gromacs, charmm]	no default
small_molecule_ff	string	Force field for small molecules. Options: espaloma, gaff, smirnoff	espaloma
protein_path	string	Path to the protein structure (PDB or PDBx). Set to null when using box_size.	null
clean_protein	boolean	Flag indicating if cleaning the protein with PDBFixer	true
keep_heterogens	boolean	Flag indicating if keeping the heterogen atoms while cleaning the protein. Waters will be always kept.	false
variants	dictionary	Dictionary of residues for which a variant is requested (different protonation state) in the form {"chain_id:res_id":"protonation_state"}, None for the rest of the residues.	{}
repulsive_forces	dict	Dict of pairwise repulsive forces: {"name": {"residueA": "BEN", "residueB": "PRP", "epsilon": 0.01, "sigma": 4.0}}. epsilon in kcal/mol, sigma in Angstrom.	{}	✅	❌	❌	❌
solvent_smiles	string	SMILES string of the solvent to use.	H2O
solvent_copies	integer	If specified, the box won't be filled up with solvent, but will have the exact number of solvent molecules specified.	null
positive_ion	string	Ion type for positive charge neutralisation.	Na+
negative_ion	string	Ion type for negative charge neutralisation.	Cl-
padding	float	Padding around the protein in Angstrom.	10.0
box_size	float	Box edge length in Angstrom. Required when no protein_path is given.	null
ligands	dict	Dict mapping ligand name to path of SDF/MOL2 file.	{}
membrane	boolean	Flag indicating if the system has a membrane or not.	false
lipid_type	string	If membrane is true specify the lipid to use. Supported lipids: [POPC, POPE, DLPC, DLPE, DMPC, DOPC, DPPC]	POPC
lipid_patch_path	string	Path to a pre-equilibrated lipid patch (mutually exclusive with lipid_type).	null
memb_cosolv_placement	string	Which side of the membrane to place the cosolvents. Options: both, outside, inside	both
waters_to_keep	list	List of residue indices of waters to preserve in membrane systems.	[]
output_dir	string	Path to the output directory for results.	null

1. Preparation

$ create_cosolvent_system -c config.yaml

2. Run MD simulations (external) CosolvKit no longer runs MD itself — that step is delegated to an external MD engine. Step 1 writes the prepared system in the format you selected (md_engine): OpenMM (system.xml + system.pdb), or Amber / GROMACS / CHARMM topology and coordinate files. Run the simulation with the engine of your choice using those files, then bring the resulting trajectory back for analysis. See examples/02_equilibration-and-production for a worked equilibration + production setup using AutoPath.

3. Analysis See examples/03_analysis for a YAML-driven analysis run (analyze_cosolvent_simulation -cfg analysis.yaml) and the interactive hotspot dashboard (visualize_hotspots).

3.1 centering, imaging, and aligning a trajectory

To generate meaningful cosolvent densities for visualization, the trajectory must be centered and aligned on the region of interest. Centering is placing a set of atoms at the center of the simulation box (without rotating the box), imaging is placing all atoms inside the box if they traveled beyond the periodic boundary conditions, and aligning is rotating and translating the system so that a selection of atoms overlaps with some reference positions.

Usually, trajectories are aligned on a macromolecule such as a protein, but parts of macromolecules that are flexible and move during the simulation can still cause densities to smear. If such flexible parts are of interest, it is a good idea to align the trajectories on each flexible part independently. If the region of interest is a specific location of a large or flexible protein, it is best to align using the vicinity of the region of interest, rather than the whole protein.

One option to align and image trajectories is cpptraj. It should be installed automatically by the installation instructions above. First we create an input file for cpptraj, which we will call process.cpptraj:

parm system.pdb
trajin trajectory.dcd
center :@CA
autoimage
reference system.pdb [myref]
rms ref [myref] @CA out protein.rmsd
trajout trajectory_aligned.xtc

There are two important selections in this input file that are system specific and need to be edited manually, the one for centering the trajectory after center command, and the one for aligning after rms. See the documentation for defining selections. To run it, system.pdb needs to be on the working directory:

cpptraj -i process.cpptraj

It will write trajectory_aligned.xtc. This trajectory should inspected to make sure the region of interest is not moving or wrapping around the periodic boundaries. First, load system.pdb into Pymol, and then type the following into Pymol's command line: load_traj trajectory_aligned.xtc, system.

An example of another program that can image and center trajectories is MDAnalysis. For imaging, see its documentation about wrapping and unwrapping.

3.2 the actual analysis

from cosolvkit.analysis import Report
"""
"""
    report = Report(statistics_file, traj_file_aligned, top_file, 
                    cosolvent_names=['BEN'], out_path=out_path)
    
    report.generate_report(equilibration=True, 
                           rdf=True,
                           rmsf=True, 
                           avg_selection='protein',
                           align_selection='protein and name CA',
                           )
    report.generate_density_maps(
                                 cosolvent_names=['BEN'],
                                 gridsize=0.375,
                                 temperature=300,
                                 use_atomtypes=True, 
                                 atomtypes_definitions='../cosolvkit/data/dacar_atomtypes.json', 
                                 )    
    report.generate_pymol_session(density_files=out_path)

Add centroid-repulsive potential

To overcome aggregation of small hydrophobic molecules at high concentration (1 M), a repulsive interaction energy between fragments can be added, insuring a faster sampling. This repulsive potential is applied only to the selected fragments, without perturbing the interactions between fragments and the protein. The repulsive potential is implemented by adding a virtual site (massless particle) at the geometric center of each fragment, and the energy is described using a Lennard-Jones potential (epsilon = -0.01 kcal/mol and sigma = 12 Angstrom).

Luckily for us, OpenMM is flexible enough to make the addition of this repulsive potential between fragments effortless (for you). The addition of centroids in fragments and the repulsive potential to the System holds in one line using the add_repulsive_centroid_force function. Thus making the integration very easy in existing OpenMM protocols. In this example, we are adding repulsive forces between BEN and PRP molecules.

from cosolvkit.cosolvent_system import CosolventSystem

cosolv = CosolventSystem(cosolvents, forcefields, simulation_format, receptor_path, radius=radius)
# build the system in water
cosolv.build(neutralize=True)
cosolv.add_repulsive_forces({"BEN_PRP": {"residueA": "BEN", "residueB": "PRP", "epsilon": 0.01, "sigma": 4.0}})

List of cosolvent molecules

Non-exhaustive list of suggested cosolvents (molecule_name, SMILES string and resname):

• Benzene C1=CC=CC=C1 BEN
• Methanol CO MEH
• Propane CCC PRP
• Imidazole C1=CN=CN1 IMI
• Acetamide CC(=O)NC ACM
• Methylammonium C[NH3+] MAM
• Acetate CC(=O)[O-] ACT
• Formamide C(=O)N FOM
• Acetaldehyde CC=O ACD

Config files

An example configuration file can be found in the data folder:

• config.yaml