Beyond Conventional MD: Exploring Rare Molecular Events

Molecular dynamics promises something close to a movie of molecular life. Atoms move, bonds flex, side chains rotate, proteins breathe, and ligands explore their molecular environment.

The problem is that many of the events we care most about in drug discovery are rare.

A conventional MD simulation can spend its entire trajectory exploring one metastable basin while a large conformational transition, ligand-binding event, or ligand dissociation remains hidden behind a free-energy barrier. The event is physically possible, but its probability within the accessible simulation time may simply be too low.

This creates a fundamental paradox: we use molecular dynamics to understand molecular motion, yet some of the most biologically important motions may never occur during the simulations we can practically run.

Over several decades, researchers have attacked this problem from different directions. Targeted molecular dynamics (TMD) drives a system toward a known structure. Random Acceleration Molecular Dynamics (RAMD) encourages a buried ligand to reveal possible exit routes. Umbrella sampling and metadynamics bias chosen coordinates to reconstruct free-energy landscapes. Replica exchange and accelerated MD, including Gaussian accelerated MD (GaMD), broaden conformational exploration. Parallel Cascade Selection Molecular Dynamics (PaCS-MD) repeatedly continues trajectories that make progress. Weighted ensemble, milestoning, and Markov state models provide routes toward statistically interpretable kinetics, while τRAMD offers an efficient way to rank relative residence times.

These methods belong to the same broad rare-event landscape, but they are not interchangeable. Each changes a different part of the computational experiment: the forces, the energy landscape, the allocation of trajectories, or the statistical analysis.

The question gradually changed from:

Can we make the event happen?

into:

Can we discover how it happens?

and eventually:

Can the sampled paths tell us something quantitative about thermodynamics or kinetics?

1994: If the endpoint is known, drive the system toward it

In 1994, Schlitter, Engels, and Krüger introduced Targeted Molecular Dynamics, or TMD.

Their starting observation was straightforward: biologically relevant conformational transitions may occur only accidentally during ordinary MD, or may require simulation conditions far from those normally used experimentally. Rather than waiting for a spontaneous transition, TMD introduces a time-dependent structural restraint that progressively drives the system toward a known target conformation.

If an open and closed experimental structure are available, for example, TMD can generate a possible pathway connecting them.

This answers an important question:

What structural sequence could connect state A to state B?

But TMD also introduced a distinction that remains central to enhanced sampling today. Because the target structure is specified in advance and the trajectory is actively driven toward it, the resulting pathway is not necessarily the pathway that would occur most frequently under unbiased equilibrium conditions.

TMD is therefore useful for generating mechanistic hypotheses, identifying moving protein regions, and exploring possible intermediates. It is not, by itself, a direct measurement of natural transition kinetics.

That distinction between finding a pathway and measuring its probability or timescale becomes increasingly important in the methods that followed.

2000: Susanna Lüdemann and the search for invisible P450 tunnels

Six years later, Susanna K. Lüdemann, Valère Lounnas, and Rebecca C. Wade confronted a different rare-event problem.

Their system was cytochrome P450cam.

P450 enzymes contain a deeply buried catalytic heme, creating an apparently simple but mechanistically important question:

How does a substrate actually reach the active site, and how does the product get back out?

Crystal structures provide static snapshots of the protein, but they do not necessarily reveal transient channels that appear only as side chains move or protein regions fluctuate. For P450cam, conventional structures did not reveal one obvious route that could explain substrate access and product egress, and spontaneous entry or exit was not something that could simply be expected to occur during accessible conventional MD.

Lüdemann and colleagues therefore introduced Random Expulsion Molecular Dynamics.

At the time, they abbreviated it REMD.

Instead of defining a destination or predetermined tunnel, the method applies an additional force to the ligand in a randomly chosen direction. If the ligand makes sufficient progress, the direction is maintained. If it becomes trapped, a new random direction is selected.

The idea is simple:

Do not tell the ligand which tunnel to use. Encourage it to escape and see which routes the protein allows.

This makes Random Expulsion MD conceptually different from TMD. TMD begins with a known structural destination. Random Expulsion MD begins with a buried ligand and asks the protein to reveal possible exits.

What the original Random Expulsion MD found

The simulations revealed three major ligand-exit pathway classes in P450cam. Importantly, one corresponded to a channel that had already been suggested from crystallographic and site-directed mutagenesis evidence, providing independent support that the accelerated simulations were identifying biologically meaningful routes.

The study also challenged another hypothesis. A water-filled channel had previously been suggested as a possible product-exit route, but ligand escape through that pathway was not observed under the simulated conditions.

The simulations further suggested that P450cam did not require a dramatic global opening of the protein to release camphor. Relatively modest backbone rearrangements, together with side-chain rotations and perturbation of local interactions, could create sufficient space for ligand passage.

That result remains conceptually important. A buried active site does not necessarily require the whole protein to swing open like a door. Transient local fluctuations can create pathways that barely exist in a static crystal structure.

Random Expulsion MD therefore helped establish an idea that is now central to molecular simulation: functional pathways can emerge from ensembles of transient protein conformations rather than from one obvious static tunnel.

A historical naming trap: REMD then is not REMD now

There is an unfortunate terminology problem for anyone reading this literature today.

In the Lüdemann paper, REMD means Random Expulsion Molecular Dynamics.

Today, however, REMD almost always means Replica Exchange Molecular Dynamics.

These are completely different methods.

Replica Exchange Molecular Dynamics runs multiple replicas of a molecular system under different conditions, most commonly different temperatures, and periodically attempts exchanges between replicas using an appropriate statistical acceptance criterion. High-temperature replicas cross barriers more readily, while exchanges help lower-temperature replicas explore configurations that would otherwise be difficult to reach. The widely cited molecular-dynamics formulation by Sugita and Okamoto appeared in 1999.

To avoid the acronym collision, Random Expulsion Molecular Dynamics became more commonly referred to as Random Acceleration Molecular Dynamics (RAMD). Modern RAMD therefore traces directly back to the method introduced by Lüdemann, Lounnas, and Wade.

When reading older literature, the context matters:

REMD in an early P450 ligand-egress paper may mean Random Expulsion MD. REMD in most modern enhanced-sampling papers means Replica Exchange MD.

Same acronym, fundamentally different purpose and sampling strategy.

From tunnel discovery to a general ligand-egress method

The significance of the Lüdemann study extends beyond P450cam.

RAMD introduced a computational strategy for a general class of problems:

• a ligand is buried inside a protein;
• the exit pathway is unknown;
• spontaneous dissociation is too slow to observe;
• one wants to discover possible routes without prescribing a specific tunnel.

Repeated trajectories can reveal recurring channels, gating residues, transient interactions, steric bottlenecks, hydration changes, and alternative molecular mechanisms.

The method therefore turns ligand dissociation from a rare accident into an event that can be sampled repeatedly.

But there is an important limitation: the extra force changes the dynamics. The nanoseconds required for a ligand to escape during RAMD are not automatically the ligand’s physical residence time.

At this stage, RAMD was primarily a pathway-discovery method.

Meanwhile, another branch of enhanced sampling focused less on pushing one ligand along a path and more on changing the landscape that makes rare events rare.

2002: Metadynamics attacks the free-energy barrier

In 2002, Alessandro Laio and Michele Parrinello introduced metadynamics, one of the most influential enhanced-sampling methods now used in molecular simulation.

Metadynamics begins by describing the event through a small set of collective variables (CVs). These might be a protein-ligand distance, a coordination number, an RMSD, a torsion, or a combination of descriptors that separates relevant states.

During the simulation, a history-dependent bias is deposited in the visited regions of CV space. In the original formulation, this bias is built from small Gaussian functions. Repeatedly visited free-energy wells are progressively filled, discouraging the system from remaining trapped and promoting exploration of new regions.

The question is therefore different from RAMD:

RAMD asks where a ligand can escape if it is encouraged to move. Metadynamics asks which free-energy barriers separate states defined by chosen collective variables.

This distinction gives metadynamics a major advantage: after appropriate reweighting and convergence analysis, it can reconstruct a free-energy surface along the selected CVs. For ligand binding or dissociation, that surface can reveal bound states, intermediates, transition regions, and relative barriers rather than only a set of exit trajectories.

Its greatest strength is also its central vulnerability. If an important slow degree of freedom is absent from the CVs—such as a hidden protein gate, ligand orientation, or water rearrangement—the system can appear stuck even while the chosen coordinate is being biased. A smooth-looking free-energy surface is not proof that the relevant motions were sampled or that the CVs were sufficient.

Standard metadynamics also distorts time. The observed transition time in a biased trajectory is not automatically a physical rate. Infrequent metadynamics was developed to recover kinetics by depositing bias rarely enough that transition-state crossings remain approximately unperturbed. This can yield rate estimates under specific assumptions, but it requires careful CV selection, convergence tests, and validation of the recovered transition-time statistics.

Other major routes across the landscape

Metadynamics is not the only important method missing from a ligand-egress-centered story. Several established approaches solve adjacent problems:

These methods should not be treated as competitors in a single ranking. Umbrella sampling is natural when the coordinate is already defensible and a potential of mean force is the target. Metadynamics is attractive when transitions and free-energy basins must be explored in a low-dimensional CV space. Replica exchange helps systems move between energetic or Hamiltonian conditions. aMD and GaMD reduce barriers without requiring the slow coordinates to be specified in advance.

For highly flexible proteins such as CYP3A4, GaMD is particularly useful for discovering receptor conformational ensembles before docking or mechanistic analysis. That is a different role from RAMD: GaMD asks which protein and ligand states become accessible when the landscape is smoothed, whereas RAMD focuses on accelerated ligand escape and route discovery.

The next method in the ligand-pathway story changes neither a single coordinate nor the potential-energy surface directly. Instead, it changes which trajectories receive more computational effort.

2013: PaCS-MD replaces forcing with selection

In 2013, Ryuhei Harada and Akio Kitao introduced Parallel Cascade Selection Molecular Dynamics, or PaCS-MD.

The rare-event problem was the same, but the solution was fundamentally different.

PaCS-MD does not continuously push the system with an additional external force. Instead, it runs many short conventional MD trajectories in parallel. At the end of each cycle, structures that have progressed furthest toward the event of interest are selected. Those structures are assigned new velocities and used to initiate another generation of short trajectories.

The procedure becomes:

short MD → selection → reseeding → short MD → selection → reseeding

until the desired state is reached.

The distinction is important. In RAMD, acceleration comes from an added force acting on the ligand. In PaCS-MD, acceleration comes primarily from which trajectories are allowed to continue.

Each individual segment can therefore be an ordinary MD trajectory without a continuous steering force. However, this does not make the overall PaCS-MD ensemble equivalent to unbiased equilibrium sampling. Repeated selection changes which trajectories survive.

Why the progress coordinate matters

Harada and Kitao demonstrated PaCS-MD using folding of the ten-residue protein chignolin and the open-to-closed transition of T4 lysozyme. Repeated cycles of short simulations were sufficient to reach structures close to the desired product states.

The approach occupies an interesting middle ground between brute-force MD and continuously steered simulations.

But its main strength also creates its main vulnerability: something must determine which trajectories are considered promising.

That means PaCS-MD depends on a progress coordinate or selection criterion. If that coordinate is poorly chosen, the method may favor structures that appear geometrically closer to the target while overlooking mechanistically important alternatives.

This is a broader lesson in enhanced sampling: accelerating sampling is only part of the problem; defining meaningful progress can be equally difficult.

When PaCS-MD was turned toward ligand dissociation

Tran, Takemura, Kuwata, and Kitao subsequently applied PaCS-MD to dissociation of tri-N-acetyl-D-glucosamine from hen egg-white lysozyme.

The contrast with conventional simulation illustrates why such methods are useful. The complex remained bound throughout a 1-μs conventional MD simulation, whereas PaCS-MD generated dissociation within roughly 1–10 ns of accumulated simulation time.

The method repeatedly selected structures showing greater ligand separation and launched new short trajectories from those states. This form is often referred to as dissociation PaCS-MD, or dPaCS-MD.

Generating a pathway also need not be the final step. The resulting states can be subjected to statistical analysis, including Markov state models, to move from trajectory generation toward free-energy or kinetic interpretation.

2018: From ligand escape to residence time

Discovering how a ligand leaves a protein is useful. Medicinal chemistry, however, often wants another quantity:

How long does the ligand remain bound?

For a simple first-order dissociation process, residence time is related to the dissociation rate by

τ = 1 / k_off

Two ligands can have similar equilibrium affinity while exhibiting very different association and dissociation kinetics. Residence time can therefore provide information that is not captured by a single equilibrium binding free energy.

The problem is computational. Ligand dissociation can occur on timescales far beyond those routinely accessible to conventional atomistic MD.

Kokh and colleagues addressed this problem using τRAMD.

Rather than interpreting one accelerated escape trajectory, τRAMD performs many RAMD simulations initiated from an equilibrated protein-ligand ensemble and analyzes the distribution of accelerated dissociation times.

The conceptual shift is from:

Which route does the ligand take?

to:

Across many accelerated exits, does one ligand consistently require longer to dissociate than another?

What τRAMD demonstrated

The original 2018 study investigated 70 drug-like inhibitors of the N-terminal domain of HSP90α.

For 78% of the compounds, the computed relative residence times were within approximately a factor of 2.3 of experiment, with still better agreement within related congeneric series.

This matters because τRAMD does not require direct simulation of the experimental dissociation timescale. Instead, it asks whether accelerated escape statistics preserve enough information about the underlying molecular barriers to rank compounds comparatively.

The trajectories can also be mechanistically informative. Transient polar contacts, steric bottlenecks, protein rearrangements, and intermediate interactions can help explain why one ligand repeatedly takes longer to escape than another.

RAMD had therefore evolved from a tool for finding hidden P450 tunnels into a method capable of addressing a question directly relevant to drug optimization.

Relative kinetics, not a magic stopwatch

The artificial force still accelerates dissociation.

A ligand leaving after 2 ns in a τRAMD trajectory does not mean its experimental residence time is 2 ns.

Instead, accelerated dissociation statistics are compared across compounds under a consistent protocol. τRAMD is therefore most naturally treated as a relative residence-time ranking method, particularly when comparing structurally related ligands.

Its reliability also depends on the starting complex, the conformational ensemble, force-field quality, RAMD parameters, and adequate replicate sampling.

This distinction is critical because enhanced sampling can easily produce visually convincing molecular movies. A convincing movie is not automatically a quantitative kinetic measurement.

When physical kinetics is the actual target

Some methods are designed from the outset to construct statistically meaningful ensembles of rare pathways rather than to generate one accelerated movie.

Weighted ensemble (WE) divides progress space into regions and periodically replicates trajectories entering undersampled regions while merging others. Their statistical weights are adjusted so that the underlying dynamics remain unbiased. With adequate sampling and correct bookkeeping, WE can estimate pathway probabilities, steady-state fluxes, and rate constants without adding a force to the system.

Milestoning divides a long transition into crossings between interfaces called milestones. Many short trajectories provide local transition probabilities and waiting times, which are then assembled to estimate quantities such as mean first-passage times. Its reliability depends on the milestone construction and on whether the assumptions used to connect local transitions are satisfied.

Markov state models (MSMs) discretize conformational space into states and estimate transition probabilities at a chosen lag time from many trajectories. They can connect distributed sampling to long-timescale populations, pathways, and kinetics, but the state definition, lag time, sampling coverage, and Markovianity must be tested rather than assumed.

Transition path sampling provides another route by sampling ensembles of reactive trajectories between states without requiring a detailed reaction coordinate in advance. It is especially valuable for mechanism discovery, although obtaining absolute rates requires additional treatment.

These frameworks are more statistically ambitious than simply observing an escape. They also demand more: clearly defined states, adequate transition sampling, convergence diagnostics, and careful uncertainty analysis.

Choose the method from the observable

This map is deliberately selective rather than exhaustive. Adaptive biasing force, steered MD, string methods, forward-flux sampling, and other approaches occupy neighboring parts of the same landscape. The key choice is not which method is most fashionable, but which observable the simulation must support.

Pathway discovery and kinetics are different scientific problems

Taken together, these methods reveal an important hierarchy.

1. Pathway generation: Can the rare event be made observable?
2. Mechanistic reproducibility: Do the same gates, contacts, or intermediates recur across trajectories?
3. Relative ranking: Can compounds be compared consistently under the same accelerated protocol?
4. Physical kinetics: Can transition statistics be connected to experimental rates and timescales?

These levels require progressively stronger assumptions and statistical treatment.

A single accelerated trajectory may reveal a plausible pathway. An ensemble can reveal recurring mechanisms. A calibrated protocol may rank compounds. But physical rate constants generally require methods specifically designed to connect transition statistics with physical time, such as Markov state models, milestoning, weighted ensemble methods, or related kinetic frameworks.

The fundamental lesson is simple:

The statistical meaning of a molecular trajectory depends on how that trajectory was generated.

Why this story is particularly relevant to CYP enzymes

Cytochrome P450 enzymes remain an especially interesting application because ligand behavior involves much more than binding affinity.

For metabolism to occur, a substrate must first reach the active site. It must sample orientations compatible with productive chemistry. A metabolically susceptible atom must approach the reactive heme species with appropriate geometry. The enzyme can reorganize around the substrate. After chemistry occurs, the product must eventually leave.

For highly flexible P450s such as CYP3A4, this produces an interconnected network of conformational dynamics, ligand access, productive orientation, chemistry, and egress.

Static docking captures only a subset of that process.

Methods such as RAMD or dPaCS-MD can provide complementary information about access channels, exit pathways, gating residues, hydration patterns, and ligand-dependent protein rearrangements.

But one distinction must remain clear:

Residence time is not metabolic turnover.

A ligand that remains bound for a long time is not necessarily metabolized rapidly. Conversely, rapid ligand escape does not automatically imply poor metabolism.

P450 catalysis additionally depends on productive orientation, accessibility of candidate sites of metabolism, intrinsic chemical reactivity, the state of the catalytic heme, and progression through the catalytic cycle.

Access and egress dynamics are therefore pieces of the metabolic mechanism, not substitutes for chemistry.

From trajectories to molecular descriptors

One particularly interesting future direction is to stop treating enhanced-sampling simulations merely as movies.

A ligand-egress trajectory contains much more information than its final escape time:

• transient protein-ligand contacts;
• gate-opening events;
• water rearrangements;
• changes in ligand orientation;
• temporary metastable states;
• residue movements;
• pathway identity;
• steric bottlenecks.

These properties can be transformed into quantitative molecular representations. Protein-ligand interaction fingerprints, pathway fingerprints, contact graphs, water-network descriptors, residue mobility features, and learned trajectory embeddings could potentially distinguish mechanisms or help predict kinetic behavior.

Machine-learning analyses of τRAMD trajectories have already illustrated this direction by using transient interaction information to help explain and predict residence-time behavior.

This creates an interesting convergence between molecular simulation and machine learning. Rather than asking ML to replace physics-based simulation, simulation can generate mechanistically informed features that are absent from conventional molecular representations.

The challenge is validation. Trajectory-derived features can be high-dimensional and protocol-dependent. Replicates, parameter-sensitivity analysis, independent test systems, and experimental kinetic data remain essential.

From accidental events to matched computational experiments

Looking back, the progression is remarkably clear.

TMD asked: How might two known protein structures be connected?

Random Expulsion MD asked: How might a ligand escape when the route is unknown?

Metadynamics asked: Which free-energy barriers separate states in a chosen collective-variable space?

aMD and GaMD asked: Can the energy landscape be smoothed without defining those coordinates in advance?

PaCS-MD asked: Can we repeatedly select promising short trajectories instead of waiting for one lucky long simulation?

τRAMD asked: Can ensembles of accelerated escapes preserve enough information to compare ligand residence times?

Weighted ensemble, milestoning, and MSMs asked: Can many short trajectories be assembled into statistically interpretable long-timescale behavior?

The methods are different, but they all attack the same underlying problem.

Rare molecular events are often too important to ignore and too slow to wait for.

The solution has therefore been to redesign the computational experiment.

Final perspective

“Accidental transitions” is a useful description of a major limitation of conventional molecular dynamics.

A protein may be capable of opening. A ligand may be capable of leaving. A substrate may have several possible tunnels. But during the finite trajectory we happen to run, none of those events may occur.

The methods discussed here replace passive waiting with different forms of designed computational exploration.

TMD guides the system toward a known endpoint.

Random Expulsion MD, now generally called RAMD, searches for unknown ligand pathways using randomly directed acceleration.

Umbrella sampling and metadynamics bias selected coordinates to reveal free-energy profiles or surfaces.

Replica exchange, aMD, and GaMD broaden exploration by changing how the system crosses energetic barriers.

PaCS-MD reallocates simulation effort toward trajectories already progressing toward a rare event.

τRAMD uses ensembles of accelerated exits to extract comparative kinetic information.

Weighted ensemble, milestoning, and MSMs use ensembles of shorter trajectories to approach physical pathways and kinetics with explicit statistical models.

Their evolution also reflects a broader maturation of molecular simulation. The goal is no longer simply to generate an impressive trajectory. The goal is to understand what the trajectory represents statistically, whether its mechanism is reproducible, and which experimental observable it can legitimately inform.

For drug discovery, that leads to three complementary questions:

Thermodynamics: How favorable is the bound state?

Mechanism: How does the system move between states?

Kinetics: How quickly do those transitions occur?

The most useful molecular-simulation workflows increasingly try to connect all three while keeping the assumptions and uncertainty of each layer explicit. There is no universally best enhanced-sampling method—only a method whose biases, statistics, and output are appropriate for the scientific question.

References

1. Schlitter, J.; Engels, M.; Krüger, P. “Targeted Molecular Dynamics: A New Approach for Searching Pathways of Conformational Transitions.” Journal of Molecular Graphics 1994, 12, 84-89. https://doi.org/10.1016/0263-7855(94)80072-3

2. Lüdemann, S. K.; Lounnas, V.; Wade, R. C. “How Do Substrates Enter and Products Exit the Buried Active Site of Cytochrome P450cam? 1. Random Expulsion Molecular Dynamics Investigation of Ligand Access Channels and Mechanisms.” Journal of Molecular Biology 2000, 303, 797-811. https://doi.org/10.1006/jmbi.2000.4154

3. Harada, R.; Kitao, A. “Parallel Cascade Selection Molecular Dynamics (PaCS-MD) to Generate Conformational Transition Pathway.” The Journal of Chemical Physics 2013, 139, 035103. https://doi.org/10.1063/1.4813023

4. Kokh, D. B.; et al. “Estimation of Drug-Target Residence Times by τ-Random Acceleration Molecular Dynamics Simulations.” Journal of Chemical Theory and Computation 2018, 14, 3859-3869. https://doi.org/10.1021/acs.jctc.8b00230

Additional papers discussed

5. Tran, D. P.; Takemura, K.; Kuwata, K.; Kitao, A. “Protein-Ligand Dissociation Simulated by Parallel Cascade Selection Molecular Dynamics.” Journal of Chemical Theory and Computation 2018, 14, 404-417. https://doi.org/10.1021/acs.jctc.7b00504

6. Sugita, Y.; Okamoto, Y. “Replica-Exchange Molecular Dynamics Method for Protein Folding.” Chemical Physics Letters 1999, 314, 141-151. https://doi.org/10.1016/S0009-2614(99)01123-9

7. Torrie, G. M.; Valleau, J. P. “Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling.” Journal of Computational Physics 1977, 23, 187-199. https://doi.org/10.1016/0021-9991(77)90121-8

8. Laio, A.; Parrinello, M. “Escaping Free-Energy Minima.” Proceedings of the National Academy of Sciences 2002, 99, 12562-12566. https://doi.org/10.1073/pnas.202427399

9. Hamelberg, D.; Mongan, J.; McCammon, J. A. “Accelerated Molecular Dynamics: A Promising and Efficient Simulation Method for Biomolecules.” The Journal of Chemical Physics 2004, 120, 11919-11929. https://doi.org/10.1063/1.1755656

10. Miao, Y.; Feher, V. A.; McCammon, J. A. “Gaussian Accelerated Molecular Dynamics: Unconstrained Enhanced Sampling and Free Energy Calculation.” Journal of Chemical Theory and Computation 2015, 11, 3584-3595. https://doi.org/10.1021/acs.jctc.5b00436

11. Tiwary, P.; Parrinello, M. “From Metadynamics to Dynamics.” Physical Review Letters 2013, 111, 230602. https://doi.org/10.1103/PhysRevLett.111.230602

12. Huber, G. A.; Kim, S. “Weighted-Ensemble Brownian Dynamics Simulations for Protein Association Reactions.” Biophysical Journal 1996, 70, 97-110. https://doi.org/10.1016/S0006-3495(96)79552-8

13. Dellago, C.; Bolhuis, P. G.; Csajka, F. S.; Chandler, D. “Transition Path Sampling and the Calculation of Rate Constants.” The Journal of Chemical Physics 1998, 108, 1964-1977. https://doi.org/10.1063/1.475562

14. Faradjian, A. K.; Elber, R. “Computing Time Scales from Reaction Coordinates by Milestoning.” The Journal of Chemical Physics 2004, 120, 10880-10889. https://doi.org/10.1063/1.1738640

15. Prinz, J.-H.; et al. “Markov Models of Molecular Kinetics: Generation and Validation.” The Journal of Chemical Physics 2011, 134, 174105. https://doi.org/10.1063/1.3565032