How to Run DomHMM

This section explains how to configure and run DomHMM with several different options.

Note

In the project’s /example directory, you can find real-life examples of DomHMM.

Running DomHMM

DomHMM’s main class is PropertyCalculation. In a basic example, it is initialized as

model = domhmm.PropertyCalculation(universe_or_atomgroup=universe,
                                   leaflet_kwargs=leaflet_kwargs,
                                   membrane_select=membrane_select,
                                   leaflet_select="auto",
                                   heads=heads,
                                   sterol_heads=sterol_heads,
                                   sterol_tails=sterol_tails,
                                   tails=tails)

Then it can be run as

model.run(start=start_frame, end=end_frame, step=step)

Main Parameters

Let’s dive into each parameter’s details.

  • In initialization process, universe_or_atomgroup parameter stands for MDAnalysis universe. It contains your simulation’s trajectory and topology file. It can be created as

path2xtc = "YOUR_XTC_FILE.xtc"
path2tpr = "YOUR_TPR_FILE.tpr"
universe = mda.Universe(path2tpr, path2xtc)

  • leaflet_kwargs parameter stands for MDAnalysis LeafletFinder function’s arguments. It is used to determine each leaflet’s lipids. leaflet_kwargs requires head groups of lipids but not sterols.

# An example where all lipids head group is PO4
leaflet_kwargs={"select": "name PO4", "pbc": True}

  • membrane_select argument is for atom group selection of universe. It is useful for simulations that contain non-membrane residues/molecules inside. If the universe contains only membrane elements, the parameter can be left in the default option which is all

# An example where simulation contains DPPC and DIPC lipids, and CHOL sterol
membrane_select = "resname DPPC DIPC CHOL"

  • leaflet_select argument is a selection option for lipids which can be a list of atom groups, a list of string queries, or automatic via LeafletFinder.

# List of atom groups
lower_leaflet = universe.select_atoms("lower leaflet lipids query")
upper_leaflet = universe.select_atoms("upper leaflet lipids query")
leaflet_select = [lower_leaflet, upper_leaflet]
# List of query strings
leaflet_select = ["lower leaflet lipids query", "upper leaflet lipids query"]
# Leave leaflet detection to DomHMM via LeafletFinder
leaflet_select = "auto"

  • heads parameter requires lipids head groups. For atomistic simulations, the head molecules’ center atom can be entered.

heads = {"DPPC": "PO4", "DIPC": "PO4"}

  • sterol_heads parameter requires sterol head groups. For atomistic simulations, the head molecules’ center atom can be entered.

# Martini Cholesterol example
sterol_heads = {"CHOL": "ROH"}
# Atomistic Cholesterol example
sterol_heads = {"CHL1": "O3"}

  • sterol_tails parameter requires sterol tail groups. It should be considered that each tail should be entered in the same order for each lipid.

# Martini Cholesterol example while ROH head as the first element and C1 start of the tail the second element
sterol_tails = {"CHOL": ["ROH", "C1"]}
# Atomistic Cholesterol example while O3 head as first element and C20 start of tail as second element
sterol_tails = {"CHL1": ["O3", "C20"]}

  • tails parameter requires lipids tail groups. It should be considered that each tail should be entered in the same order for each lipid.

# Example of tails in order of {"Lipid_1":[[Acyl_Chain_1],[Acyl_Chain_2]], "Lipid_2":[[Acyl_Chain_1],[Acyl_Chain_2]]}
tails = {"DPPC": [["C1B", "C2B", "C3B", "C4B"], ["C1A", "C2A", "C3A", "C4A"]],
             "DIPC": [["C1B", "D2B", "D3B", "C4B"], ["C1A", "D2A", "D3A", "C4A"]]}

  • For run option, you can have start, stop and step options. These options arrange which frame to start or stop. You can also set the model to be trained for each X frame by setting step=X.

# An example where DomHMM model training starts from the 5th frame and ends in the 1000th frame while taking each 5th step. The first three frames will be the 5th, 10th, and 15th frames.
model.run(start=5, stop=1000, step=5)

Warning

If detailed post-analysis will be conducted on results such as usage of Getis_Ord results, the input order of lipids and sterols should be in same order as in simulation. If simulation lipids are in order of DPPC, DIPC, CHOL with respect to residue IDs, keys of heads, tails, sterol_heads, and sterol_tails should be in the same order just like in this example.

Note

Since DomHMM uses the Gaussian Mixture Model and Gaussian-based Hidden Markov Model, it is suggested to not use too short or too long simulations. Short simulations may not create sensible results and long ones would take too much time to train the model. In our examples, we used simulations that contains around 2000 frames, and the model run is finished around 25-30 minutes.

Optional Parameters

  • asymmetric_membrane

It needs to be enabled if leaflets are not symmetric. With this option, models are fitted by separated data for each leaflet.

  • do_clustering

Whether to perform the hierarchical clustering or not (Default is True).

  • result_plot

Plotting option for debugging. While enabled, DomHMM will print Hidden Markov model convergence, prediction results, Getis-Ord statistic results, and clustering results of three frames.

  • save_plots

Option for saving result plots in pdf format.

  • verbose

Verbose option for debugging. It shows which steps are done in the analysis.

  • lipid_leaflet_rate

The frame rate for checking lipids leaflet assignments via LeafletFinder. In the default option, it is equal to 0 which means leaflet assignment is only done at the beginning of the analysis.

  • sterol_leaflet_rate

The frame rate for checking sterols leaflet assignments via LeafletFinder. In the default option, it is equal to 1 which means sterols leaflet assignment will be calculated in every time frame to capture flip-flops.

  • gmm_kwargs

Parameter option for Gaussian Mixture Model training. An example of this is

gmm_kwargs = {"tol": 1E-4, "init_params": 'random_from_data', "verbose": 0,
                  "max_iter": 10000, "n_init": 20,
                  "warm_start": False, "covariance_type": "full"}

  • hmm_kwargs

Parameter option for Gaussian-based Hidden Markov Model training. An example of this is

hmm_kwargs = {"verbose": False, "tol": 1E-4, "n_iter": 1000,
                  "algorithm": "viterbi", "covariance_type": "full",
                  "init_params": "st", "params": "stmc"}

  • trained_hmms

Parameter option for reusing past DomHMM HMM models. If there are several analyses that will be conducted with slightly different membrane simulations or with different parameter options, the first analysis HMM model can be reusable with this parameter.

model.run()
with open(f'hmm_model_dump.pickle', 'wb') as file:
    pickle.dump(model.results["HMM"], file)
...
with open(f'hmm_model_dump.pickle', 'rb') as file:
    reuse_hmm_model = pickle.load(file)
model_2 = domhmm.PropertyCalculation( ... ,
                                     trained_hmms=reuse_hmm_models)

  • n_init_hmm

Number of repeats for HMM model training. HMM models can be trained multiple times to achieve better performance.

  • frac

The fraction of box length in x and y outside the unit cell is considered for area per lipid calculation by Voronoi. It is an optimization process parameter that is set to 0.5 as the default.

  • p_value

Probability value that is used for z-score calculation. It is a determination percentage for domain identification with the Getis-Ord statistic. In default, it is set to 0.05 or %5.

  • tmd_protein_list

Transmembrane domain (TMD) protein list to include area per lipid calculation. TMD proteins take up space in the exoplasmic, cytoplasmic leaflets. Three backbone atoms of protein that are in close position to lipid head groups should be included in this parameter to increase the success of identification.

# Selecting three backbone atoms that are touching the exoplasmic leaflet
upBB = uni.select_atoms('name BB')[0:3]
# Selecting three backbone atoms that are touching the endoplasmic leaflet
loBB = uni.select_atoms('name BB')[-3:]
# List can be expended with multiple dictionary objects as in more than one TMD protein scenario.
tmd_protein_list = [{"0": upBB, "1": loBB}]

We encourage you to check Tips for Usage section which may contain useful information for your progress.