What if molecular dynamics could be treated like language: structures become tokens, trajectories become sentences, and every generated path comes with an explicit probability? This paper tests that idea across proteins, temperatures, and force-field conditions.

Protein conformational populations and rates of interconversion depend on physical conditions and underlying molecular interactions. Molecular dynamics is computationally demanding because atomistic simulations use femtosecond integration timesteps, while conformational transitions may take microseconds to milliseconds or longer.

Enhanced sampling methods such as metadynamics and Hamiltonian replica exchange accelerate rare-state sampling, but remain costly and system-specific because they rely on predefined collective variables or simulated thermodynamic states. Data-driven models offer an alternative by learning statistical descriptions of molecular ensembles and dynamics from existing data.

Markov state models and VAMPnets infer conformational kinetics, while BioEmu generates equilibrium-like conformational ensembles without explicit time evolution. TITO and PLaTITO learn finite-lag transition laws and have been shown to recover equilibrium probabilities or folding kinetics, but their continuous-model transition likelihoods typically require probability-flow ODE integration at each step.

MD-LLM treats molecular dynamics as a language-modeling problem. Structural tokens play the role of words, conformations the role of sentences, and trajectories the role of documents; the model learns how one conformation follows another from molecular history. Because every token comes from a normalized categorical distribution, the probability of a transition and of a complete tokenized trajectory follows directly from multiplying token probabilities, without numerical integration.

MD-LLM-2 extends the protein-specific approach to a transferable model trained across protein domains, while retaining frame-by-frame autoregressive generation with ESM3 structural tokens. The model conditions its finite-lag transition distribution on molecular history, sequence, physical temperature, and an optional Hamiltonian coupling.

Unlike a model using only the current reduced state, MD-LLM-2 explicitly conditions on recent frames, allowing some earlier memory to be represented rather than absorbed into the lag choice. The same autoregressive formulation is applied at a zero point two picosecond lag for alanine dipeptide and a one nanosecond stored-frame lag for mdCATH.

A tokenizer maps each configuration to structural tokens, and the next token frame represents the configuration one lag later while collecting sequence, temperature, and optional coupling as conditioning variables. Given the tokenized history, the model generates the next frame one token at a time.

Each token is drawn from a normalized categorical distribution over structural codes, so a fixed starting frame and sampling policy define a normalized distribution over token paths. The likelihood of any generated or molecular-dynamics trajectory can be evaluated directly.

When the token space is small, the distribution can be enumerated and state populations can be summed exactly without sampling. The predictive distribution at each step also exposes model uncertainty, which can guide sampling. Frame-isoentropic sampling uses predictive uncertainty to limit rollout drift caused by exposure bias, where generated outputs replace training data in the autoregressive history.

In T4L rollouts, mean categorical entropy rises with structural drift, with a median within-rollout Spearman correlation of zero point eight five five. After each generated frame, a controller selects one sampling temperature for the next frame so that mean predictive entropy returns to a reference value set from probes of the first transition.

The feedback acts on token probabilities while the physical temperature and Hamiltonian coupling remain fixed. MD-LLM-2 was trained on mdCATH trajectories of protein domains up to 200 residues, simulated at five temperatures between 320 and 450 kelvin.

Each frame is represented by one ESM3 structural token per residue, and the model conditions on sequence, physical temperature, the current frame, and up to three earlier frames. Without further training, evaluation covered 12 unseen mdCATH domains, Trp-cage for temperature response, and L99A T4L for rare conformational states.

For Hamiltonian conditioning, alanine dipeptide was simulated in implicit solvent with torsional energy terms scaled by a coupling lambda; lambda equal to one recovers the original force field, while lambda equal to zero removes those terms. Training used Hamiltonian replica-exchange simulations from lambda zero point eight to one point three and fixed-lambda molecular dynamics at lambda equal to one point zero, one point three, and two point zero.

Independent references at excluded couplings came from a separate Hamiltonian replica-exchange ladder for equilibrium populations and fixed-lambda molecular dynamics for first-passage kinetics. At lambda values from zero point two to zero point five, MD-LLM-2 followed the reference population trend, with a small systematic overestimate, and was closer to the references than MBAR.

Because each frame had two tokens from 128-code vocabularies, all 128 squared code pairs could be enumerated and their basin-one probabilities summed. Across four unseen couplings, ten independently trained models reached a mean absolute error of zero point zero nine one eight, compared with zero point three four nine five for MBAR reweighting.

Every individual model improved on MBAR, with errors ranging from zero point zero three five eight to zero point three two one seven. At lambda equal to zero, frame-isoentropic sampling brought MD-LLM-2 close to the reference population. Across 200 generated frames, the basin-one population was zero point three seven seven, compared with zero point three nine six in independent molecular dynamics.

Fixed sampling at tau equal to zero point seven gave zero point five three seven, while MBAR gave zero point zero five eight. The controller set its entropy target from the model's own predictions at lambda equal to zero and used no reference data at that coupling. In panel a, the blue ensemble curve stays close to the black independent-MD points across the four unseen couplings, while the orange MBAR curve sits visibly higher, especially toward the lower coupling values.

Panel b shows substantial spread among the ten independently trained models, with several points near the ensemble line but others markedly farther away; the MBAR reference line is above every model point. At held-out couplings, MD-LLM-2 also reproduced finite-time kinetics. Using autoregressive generation at a zero point two picosecond lag, the model produced 100-picosecond paths at lambda equal to one point one five and zero point eight, and these were compared with direct molecular dynamics at 20, 50, and 100 picoseconds.

Mean absolute error across both couplings, both directions, and all three horizons was zero point zero eight one one, compared with zero point one five six eight for a reversible Markov state model. The largest difference arose for transitions into basin one, which the Markov state model predicted to occur too quickly.

In the left column, the teal MD-LLM-2 curves stay closer to the direct-MD points than the dashed reversible-MSM curves, especially at longer horizons in both coupling conditions. On 12 unseen mdCATH domains, generated trajectories reproduced the overall size of the molecular-dynamics ensemble.

For each domain, one 200-frame trajectory was generated with frame-isoentropic sampling from the first frame and compared with the original trajectory. The median C-alpha radius of gyration was within one angstrom of the molecular-dynamics median in 10 of 12 domains, with a median absolute difference of zero point three nine angstroms.

In panel a, the teal and gray traces often overlap in their compact-state fluctuations, while several domains show teal excursions that the reference traces do not mirror closely. In panel b, the points cluster near the dashed equality line and within the shaded one-angstrom band, whereas panel c places every point below equality, showing a systematic separation in native-contact Jaccard medians.

Frame-isoentropic sampling corrected drift most where drift was strongest. In a model trained without history, the 200-frame median Jaccard index increased in all ten held-out 320-kelvin trajectories, with a median gain of zero point zero nine five.

In the history-conditioned model, the effect on 12 unseen domains was concentrated in the largest departures from molecular dynamics. The mean absolute difference in median C-alpha RMSD fell from 3.6 angstroms with standard sampling to 2.3 angstroms with frame-isoentropic sampling, while the largest difference fell from 11.5 to 4.2 angstroms.

In the bottom row, the dashed standard curves fall well below the blue frame-isoentropic curves, especially after roughly frame one hundred. The insets show the clearest separation in panel d, where the displayed medians are zero point four zero five for frame-isoentropic sampling and zero point two seven five for standard sampling; panel c shows a smaller gap.

For Trp-cage, an external system not included in mdCATH, compaction decreased as the conditioned temperature rose. Five trajectories were generated at each of five training temperatures from 320 to 450 kelvin, using the same starting structure and a compact-frame criterion based on radius of gyration and C-alpha contacts.

Compact occupancy decreased monotonically from 35.4 percent at 320 kelvin to 12.0 percent at 450 kelvin. At 320 kelvin, extended 1,000-frame trajectories made repeated reversible excursions from compact conformations, with compact frames making up 20 to 35 percent of frames 501 to 1,000 in every trajectory.

Panel c makes the intermittency visible: the line repeatedly rises toward roughly ninety percent per window and falls to zero, rather than staying at a fixed compact fraction. Conditioning on recent frames kept the T4L scaffold intact during generation.

From the same starting structure and five seeds, the median C-alpha RMSD was 3.4 angstroms with earlier-frame history, compared with 12.0 angstroms when history slots used the current frame and 15.5 angstroms for a model trained without history. Frame-isoentropic sampling kept the median at 3.4 angstroms and reduced the 90th-percentile RMSD from 6.3 to 4.9 angstroms, improving four of five trajectories.

The two rows make clear that the sampling comparison affects the upper tail of the RMSD distribution, not its reported median. Within the stable T4L scaffold, generated trajectories reached the transition region repeatedly. Ten 200-frame trajectories were generated with history and frame-isoentropic sampling, and each frame was assigned to the closest published structure using locally fitted backbone RMSD.

Because the model generates backbone coordinates only, the assignment reflects backbone changes such as Phe114 psi-angle rotation and helix merging, not burial of the Phe114 sidechain itself. In seven of ten trajectories, transition-region structures became the nearest reference with local RMSDs from 0.64 to 1.03 angstroms, and trajectories returned to E-like states between visits.

Four trajectories placed the transition region within 3 angstroms after global scaffold alignment, while two trajectories reached frames nearest to the excited-state B family after first visiting the transition region. In panel b, the history-based plots show the state-8-to-9 transition region appearing intermittently across the trajectory, while the no-history model remains concentrated near a narrower set of states before shifting toward red, higher-RMSD assignments.

Panels c and d show why a frame can look globally displaced yet locally match the transition region closely: the locally fitted overlay is visibly tighter around residues 100 to 120 than the whole-protein alignment. Molecular trajectories can be represented as conditional distributions over tokenized paths, conditioned on molecular history and thermodynamic state.

Because each path has an explicit probability, likelihood ratios between models and reweighting of generated paths can be evaluated directly without estimating continuous configuration-space densities. The protein model generates backbone coordinates only, so sidechain-dependent events such as Phe114 burial can be inferred only from backbone signatures; an all-atom representation would resolve those events directly.

Hamiltonian conditioning was demonstrated only in alanine dipeptide in this work. Because mdCATH varies temperature but not the Hamiltonian, extending Hamiltonian conditioning to proteins will require datasets with controlled Hamiltonian variation.

The protein lag is fixed by the mdCATH stored-frame interval, and generated protein-trajectory timescales have not yet been calibrated against molecular dynamics. The kinetic test is limited to first-passage probabilities between two alanine-dipeptide basins over 100 picoseconds.

Quantitative assessment of rates and timescales across increasingly complex systems, including proteins, is ongoing work. MD-LLM-2 shows that a transferable language model can reproduce several structural, thermodynamic, and limited kinetic behaviors while assigning explicit probabilities to molecular paths.

The result is promising, but protein timescales and broader Hamiltonian conditioning remain open challenges.