How do I setup my MSTIS simulation correctly?

Question

I'm trying to setup a Multi State Transition Interface Sampling (MSTIS) simulation to study the ligand binding of my enzyme using OPS and GROMACS. My enzyme acts on small sugar chains and I would like to simulate which sugar unit it bound most often (roughly depicted in the picture below). I tried to run a Transition Path Sampling (TPS) simulation, but so far I was not able to obtain a decorrelated path:

According to my understanding so far, TIS might work better here since I think that this transition has a quite small probability. And if I include the other bound states as well, I'm hoping to get a good correlation between the transition frequencies here and those I observe in the lab.

I docked my substrate in two ways, once with the second (AaAA) and once with the third (AAaA) sugar unit (short A) bound to the active site (the sugar at the active site is indicated by a lower case a). Using gromacs, I created an initial trajectory for both pulling them away from the active site.

I load these trajectories as described by dwhswenson here. As the collective variables (CVs), I chose the distance of one atom in the corresponding sugar unit to the metal ion in my active center (dist_AAaA & dist_AaAA).

Below a distance of 0.25, I consider the ligand to be bound, while with a distance greater than 1.5 (for both distances), I consider the ligand to be unbound. Since the distance is almost never smaller than 0.2 and the two distances influence each other, the two states (bound_AAaA & bound_AaAA) are quite small in the plot shown above. The dotted lines should indicate the interfaces.

So far, my code looks like this:

import openpathsampling as paths
from openpathsampling.engines import gromacs as ops_gmx

import mdtraj as md
import numpy as np

# define the options for the MD simulation
options = {
    'gmx_executable': 'gmx -nobackup ',
    'snapshot_timestep': 0.02,
    'n_frames_max': 200,
    'grompp_args': '-n index.ndx',
    'mdrun_args': '-nb gpu -nt 6'
}

# define the file names
trr_file_AAaA = "md_pull_AAaA_1st.trr"
trr_file_AaAA = "md_pull_AaAA_1st.trr"
pdb_file = "npt_AAaA_conf2.pdb"
top_file = "topol_AAaA+AaAA.top"

# setup the engine
TIS_engine = ops_gmx.Engine(gro=pdb_file,
                             mdp="TPS.mdp",
                             top=top_file,
                             options=options,
                             base_dir=".",
                             prefix="TIS").named("TIS_engine")

engine_setup = TIS_engine.current_snapshot

# define the CVs for all 3 states
dist_AAaA = paths.MDTrajFunctionCV(
    "dist_AAaA", md.compute_distances, engine_setup.topology, 
    atom_pairs=[[3136,3180]], periodic=False)

dist_AaAA = paths.MDTrajFunctionCV(
    "dist_AaAA", md.compute_distances, engine_setup.topology, 
    atom_pairs=[[3136,3207]], periodic=False)

dist_unbound = paths.MDTrajFunctionCV(
    "dist_unbound", md.compute_distances, engine_setup.topology, 
    atom_pairs=[[3136,3192]], periodic=False)

# define the states
bound_AAaA = (paths.CVDefinedVolume(
    dist_AAaA, 0, 0.25)).named("bound_AAaA")
bound_AaAA = (paths.CVDefinedVolume(
    dist_AaAA, 0, 0.25)).named("bound_AaAA")
unbound = (paths.CVDefinedVolume(
    dist_unbound, 1.5, 2)).named("unbound")

# define the interfaces
bound_AAaA_interface = paths.VolumeInterfaceSet(
    dist_AAaA, 0, [0.25, 0.35, 0.45, 0.55])
bound_AaAA_interface = paths.VolumeInterfaceSet(
    dist_AaAA, 0, [0.25, 0.35, 0.45, 0.55])
unbound_interface = paths.VolumeInterfaceSet(
    dist_unbound, 2, [1.5, 1.3, 1.1, 0.9, 0.8, 0.7, 0.6])

# import the two inital trajectories
# calculate the number of frames in both input trajectories
n_frames_AAaA = len(md.load(trr_file_AAaA, top=pdb_file))
n_frames_AaAA = len(md.load(trr_file_AaAA, top=pdb_file))
# external_traj uses externally-stored snapshots
init_traj_AAaA = paths.Trajectory(
    [TIS_engine.read_frame_from_file(trr_file_AAaA, num)
     for num in range(n_frames_AAaA)])
init_traj_AaAA = paths.Trajectory(
    [TIS_engine.read_frame_from_file(trr_file_AaAA, num)
     for num in range(n_frames_AaAA)])

# create two networks to obtain the initial trajectories
tps_network_for_init_traj_AAaA = paths.TPSNetwork.from_states_all_to_all(
    [bound_AAaA, unbound])
tps_network_for_init_traj_AaAA = paths.TPSNetwork.from_states_all_to_all(
    [bound_AaAA, unbound])

# take the subtrajectories matching the ensemble 
subtrajectories = []
for ens in tps_network_for_init_traj_AAaA.analysis_ensembles:
    subtrajectories += ens.split(init_traj_AAaA)

for ens in tps_network_for_init_traj_AaAA.analysis_ensembles:
    subtrajectories += ens.split(init_traj_AaAA)

# create a TIS network for the sampling itself
mstis_network = paths.MSTISNetwork(
    [(bound_AAaA, bound_AAaA_interface),
     (bound_AaAA, bound_AaAA_interface),
     (unbound, unbound_interface)])

# define a move scheme for the TIS simulation
scheme = paths.OneWayShootingMoveScheme(mstis_network, 
                                        selector=paths.UniformSelector(),
                                        engine=TIS_engine)

# make subtrajectories into initial conditions (trajectories become a sampleset)
initial_conditions = scheme.initial_conditions_from_trajectories(subtrajectories)

# check that initial conditions are valid and complete
# (raise AssertionError otherwise)
scheme.assert_initial_conditions(initial_conditions)

# setup the sampler
sampler = paths.PathSampling(
    storage=paths.Storage("AAaA_AaAA_TIS.nc", "w", engine_setup),
    move_scheme=scheme, sample_set=initial_conditions)

# sampler.run_until_decorrelated()
# run only 100 simulations for testing purposes
sampler.run(100)

sampler.storage.close()

The code runs without any error, but I'm not sure if the setup is correct since my resulting path tree looks like this:

As you may have noticed, I didn't setup the unbound state as described above since a certain CV is needed for the interface setup (unbound_interface = paths.VolumeInterfaceSet(CV,...)). Here I used the distance to another atom in the center of my ligand. I get how you would setup the unbound state depending on the two distances with

unbound = (paths.CVDefinedVolume(dist_AAaA, 1.5, 2) &
           paths.CVDefinedVolume(dist_AaAA, 1.5, 2)).named("unbound")

but I don't know how to incorporate that into the interface setup. Preferably I would like to not use an unbound state at all since (at least from my understanding) it's not a stable state and I'm mainly interested in the transition towards to two bound states crossing the outer interface first and finally ending in the stable bound state. But I'm not sure if you could set it up like this. Maybe someone could comment on this.

Furthermore, I'm not sure if the rest of the setup is correct. I did not define an ms_outers for the paths.MSTISNetwork as I'm not sure what this would be in my case or how to correctly define it. Maybe this could be used to remove the unbound state?

And I would like to specify that I'm only interested in the transitions from the unbound to either of the bound states but not in the transition from one bound state to another bound state. So something similar to this: tps_network = paths.TPSNetwork([unbound, bound_AAaA],[unbound, bound_AaAA])

I hope to get some help on these main issues, but if you have any comment regarding this setup or a link to a related paper this would be highly appreciated as well.

Answer 1

This is a great question, with a lot to unpack in it. (In other words, strap in for a long answer.) There's a bit of an XY problem here: the specific question asked is a good one, but I think the more fundamental issue needs a different solution. I'll start with the fundamental issue, then move on to the specific question.

When to use TPS and when to use TIS

In general, the choice of standard transition path sampling (TPS) vs. transition interface sampling (TIS) doesn't depend on the physical system you're modeling. Instead, it's about what information you want to calculate about the system. TPS mainly gives you mechanisms. TIS gives you mechanisms and reaction rates. The downside is that TIS is a lot more expensive than TPS, so you only want to use TIS if you need that more detailed information.

What to do when you're not getting decorrelated trajectories

The OP's plot of the path tree is a perfect example of a common issue with getting decorrelated trajectories with one-way shooting TPS. (In two-way shooting, the same problem manifests as shooting moves always getting rejected.)

Let me zoom in on part of the path tree. Under it, I'm going to sketch my guess of what the committor probability looks like as a function of this path.

Reminder: the committor is the probability that, from that configuration with random velocities, you end up in the final state before hitting the initial state. A shot in one-way shooting isn't exactly a committor attempt, because the velocities aren't truly randomized, but for condensed phase systems the velocity memory is relatively short, so you can think of each shot as "almost" a committor attempt.

The problem here is every trial from a shooting point left of the committor jump has near-0% probability of accepting a forward shot (you only accept backward shots), and every trial right of the committor jump has a near-0% probability accepting a backward shot. So you end up with this gap between blue (backward) and red (forward), and that gap means no decorrelated trajectories.

Now that we know the problem, the answer (in principle) is easy: select shooting points from closer to the transition state (where the committor is 0.5 and you have equal chance of accepting forward or backward shots). Unfortunately, easy in principle doesn't always mean easy in practice.

• If you know where the transition state is, then instead of using uniform shooting point selection, you can select from points more likely to be near the transition state. In OPS, you do this by changing the selector parameter, in, e.g., the OneWayShootingMoveScheme. There are two main options here that are relevant: BiasedSelector, which allows you to use any function you want to for the shooting point bias, and GaussianBiasSelector, which specifically uses a Gaussian form (since that's the most common choice).

• If you don't know where the transition state is, spring shooting will find its way to the transition state. In OPS, you use a different move scheme to do spring shooting. So instead of the standard TPS OneWayShootingMoveScheme, you use SpringShootingMoveScheme. Here's an example of spring shooting used on a toy model.

• Caveat: Multiple states make this harder. This specific system has a transition that, at some point, splits into two different channels (depending on which type of bound state). The risk is that if the transition state comes after the split into two channels, then you have different transition states for each channel. Path sampling might get stuck sampling only one channel. Approaches to enhance switching between channels is still an active area of research. (If your transition state is before the channel split, then will shoot from the part of the path that is common to both channels and you should see channel switching.)

Conceptually, a similar analysis was done in Arjun et al., although in that case the different states (two kinds of solids) were only identified after the sampling was done. However, that can give some ideas of how path density analysis can be useful here.

Selecting only certain state pairs in TPS

You have two bound states and one unbound state, and you're interested in unbound -> bound transitions. The standard constructor for TPSNetwork takes parameters initial_states and final_states, so you can get the network you want with:

network = paths.TPSNetwork(initial_states=[unbound],
                           final_states=[bound_AaAA, bound_AAaA])

The from_states_all_to_all constructor is literally just a shortcut so that, if the initial states are the same as the final states, you don't have to type the list of states twice.

Not immediately relevant to OP, but for a more complicated set of transitions where you only wanted some transition pairs selected, you could use the from_state_pairs constructor. For example, imagine you have states A, B, C, and D, and you want transitions A->C, B->C A->D, but not B->D. from_states_all_to_all gives all possible 2-state pairs from those 4 states, (except self-transitions, like A->A). The normal constructor TPSNetwork(initial_states=[A, B], final_states=[C, D]) would allow the B->D transitions. But you could get exactly the desired transitions with:

network = paths.TPSNetwork.from_state_pairs([(A, C), (B, C), (A, D)])

Where MSTIS isn't doing what you might want

There are a couple of points to make about your MSTIS setup -- first, as I mentioned earlier, unless you want something that you can calculate efficiently with TIS (e.g., rates), there's no reason to do the more expensive setup. Also, I don't think the setup/results here are really doing what you want:

1. TIS allows A->A trajectories. This is necessary for the way it calculates rates, but if you're mainly interested in the mechanism, that's pointless. The short trajectories you see are A->A trajectories.

2. MSTIS links all states to all other states: This means that with MSTIS, you are allowing the AaAA to AAaA transition (and vice versa). Multiple interface set TIS (MISTIS) is an approach (disclaimer: I developed) to avoid that problem. Here's an OPS example of MISTIS.

3. The path tree follows a single replica: The path tree shown doesn't tell us much because it is probably only following the innermost replica. In MSTIS, you'll be sampling one replica for each interface (and if that path tree is for the replica 0, it's probably an innermost interface anyway, so unlikely to see transitions).

   That said, it looks to me like this is an example of what I call the 3-frame problem in TIS -- after step 45, these trajectories are all 3 frames long, which means only one frame is not in the state (only one valid shooting point). These can be very hard to decorrelate.

On setting up the interface set/state definitions for the unbound state

First, your unbound state is a stable state -- it is stable due to entropy. There are a lot of ways for your ligand to be unbound, not so many for it to be bound. And if you just run MD in the unbound state, presumably it takes a long time before you reach the target bound states.

More generally, as discussed above, you probably don't actually want to do TIS here. But a couple points if you did:

1. Interface sets should be in terms of increasing CV values. The unbound interface set uses decreasing values of the CV. Although OPS allows you to set up interface sets where the interface values decrease, the analysis tools can't handle that (yet). The easy trick here is to wrap your distance CV in another CV that just returns the negative of it:

def neg_cv(snapshot, cv):
    return -cv(snapshot)

neg_dist_unbound = paths.CoordinateFunctionCV(
    'neg_dist_unbound',
    neg_cv,
    cv=dist_unbound
)

2. The state is treated separately from the interface set. States do not have to be defined only in terms of the CV used in the interface set. If your state definition is more strict than your innermost interface, then you have an "interstitial" region between the state and the innermost interface. On the other hand, if your state definition is less strict that your innermost interface, then the edge of your innermost interface is effectively defined by your state. All of this is fine mathematically, although it can make it harder to think about what's going on (especially if the effective interface is defined by the state definition in some degrees of freedom and the interface CV in another.)

3. You can use a custom CV that combines your existing CVs. The CV associated with the unbound interface set should be seen as "pushing you away from" the unbound state. One way to do that (independent of which final state it goes to) is to create a CV that uses both of the individual distances. For example:

def neg_min_dist(snapshot, cv1, cv2):
    # using negative here for reason listed above
    return -min(cv1(snapshot), cv2(snapshot))

unbound_cv = paths.CoordinateFunctionCV(
    "unbound_cv",
    neg_min_dist,
    cv1=dist_AaAA,
    cv2=dist_AAaA
)

Additional notes

• I notice that you're passing the periodic=False argument to MDTraj in your CVs. Since you're explicitly doing that, I assume you mean to, but I'm not sure why. It looks to me like a potential source of errors.
• Python supports float('infinity') -- for unbound states I usually use that for the maximum, because that way I can reuse the same state definition regardless of things like box size.