The emergence of the first stars marked a fundamental transition in the universe's evolution, as these stellar bodies forged the first heavy elements, altered the chemistry of primordial gas, and catalyzed the transformation of a nearly dark cosmos into the structured universe we observe today. Reconstructing this transition through numerical modeling remains a significant computational undertaking. By utilizing the MEGATRON suite of high-resolution cosmological radiation-hydrodynamics simulations, researchers are leveraging supercomputing to investigate how the chemical signatures of the universe’s first stars may be preserved within the smallest galaxies orbiting the Milky Way. 

This research underscores a critical advancement for the high-performance computing community: accurately modeling the earliest stages of galaxy formation requires simulations that simultaneously resolve gravity, hydrodynamics, adaptive mesh refinement, radiative transfer, stellar formation, stellar feedback, and non-equilibrium chemistry. MEGATRON integrates these complex physical processes at unprecedentedly small spatial scales. 

Consequently, this computational experiment tracks the evolution of a proto-Milky-Way environment from the early universe to a redshift of approximately z ≈ 8, resolving structures on a parsec scale. Ultimately, these simulations provide a robust numerical laboratory for addressing a foundational inquiry: how did the first stars dictate the chemical composition of the smallest galaxies present in the modern universe?

A universe inside the supercomputer

MEGATRON is not a conventional galaxy simulation that represents a large volume of the universe at relatively coarse resolution.

It uses a zoomed cosmological simulation focused on a region that will ultimately develop into a Milky Way-mass system.

The dark-matter particle mass is approximately 2.5 × 10⁴ solar masses, while the spatial resolution evolves from roughly 2.5 parsecs at z = 25 to approximately 5 parsecs at z = 8.

That resolution is critical.

The researchers are attempting to resolve the environments where Population III, or first-generation, stars form inside some of the earliest and smallest galaxies.

Instead of simply prescribing the effects of the first stars, the simulations can follow much of the physical environment in which those stars emerge.

The calculations use the RAMSES-RTZ adaptive-mesh-refinement code, allowing computational resolution to be concentrated where the underlying physical processes demand it.

That is an important HPC strategy.

A uniform grid capable of providing parsec-scale resolution across an entire cosmological volume would be computationally prohibitive. Adaptive mesh refinement instead allows the simulation to devote increasingly fine computational resources to the dense, dynamically important regions where stars and galaxies are forming.

The result is a multiscale problem: cosmological structure must be followed over enormous distances while the simulation simultaneously resolves processes occurring on scales measured in parsecs.

More than hydrodynamics

The computational difficulty becomes even clearer when the physical model is examined.

MEGATRON simultaneously incorporates radiative transfer and non-equilibrium chemistry involving more than 80 primordial species, molecules and metal ions.

That matters because radiation from the first stars changes the chemical and thermal state of surrounding gas.

The chemistry then affects cooling.

Cooling changes gas collapse.

Gas collapse changes star formation.

Stars generate radiation and mechanical feedback.

Supernovae distribute newly created elements through the surrounding medium.

Those metals subsequently influence later generations of stars.

In other words, the simulation is not executing one isolated physical calculation.

It is solving a tightly coupled chain of processes in which the output of one physical system becomes the input to another.

For HPC systems, that is precisely the kind of workload that makes scientific simulation difficult: the computational challenge is not simply the number of particles or grid cells, but the interaction of many physical models operating simultaneously across dramatically different spatial and temporal scales.

Four universes, one experiment

MEGATRON is actually a simulation suite rather than a single calculation.

The researchers use the same initial conditions for four simulations, varying the Population II star-formation and feedback models while keeping their Population III assumptions fixed.

The four models include different treatments of star formation and feedback, allowing the researchers to determine which features of the resulting dwarf-galaxy population remain robust when the uncertain physics of later stellar generations is changed.

That makes the project particularly interesting from an HPC perspective.

The supercomputer is being used not merely to produce one enormous dataset, but to perform a computational experiment across multiple physical models.

The researchers can then compare the resulting galaxy populations and determine which observable signatures arise repeatedly.

Across the simulations, the calculation produces a population of more than 500 simulated dwarf galaxies.

That statistical population is essential.

A single simulated galaxy could produce an interesting result. Hundreds of simulated galaxies allow researchers to ask whether a phenomenon is a systematic consequence of the underlying physics.

The iron mystery

The scientific target is an unusual feature observed in ultra-faint dwarf galaxies around the Milky Way.

At very low stellar masses, roughly below 10⁵ solar masses, observed dwarf galaxies show an approximately flat iron-metallicity relation centered around:

⟨[Fe/H]⟩ ≈ −2.5

The origin of that plateau has been debated.

Several possibilities have been proposed, including variations in Population II star formation, Population III physics, external enrichment, and the efficiency with which supernova products escape dwarf galaxies.

MEGATRON provides a computational route to test those possibilities.

The simulations naturally produce an iron-metallicity plateau at low stellar masses, broadly consistent with observations.

The researchers attribute the effect primarily to the chemistry created by massive Population III stars.

And that is where the computational experiment becomes particularly interesting.

One explosion can change an entire galaxy

The simulations indicate that most of the faint dwarf galaxies undergo approximately one Population III explosion before transitioning to Population II star formation.

The massive first-generation star can explode as a pair-instability supernova, producing substantial quantities of heavy elements.

But the supernova does not simply blow all of those elements into intergalactic space.

In sufficiently massive early halos, the gravitational potential can retain the enriched material.

The next generation of stars then forms from gas carrying the chemical signature of the first explosion.

MEGATRON finds that Population III stars tend to form in halos with masses of at least approximately 5 × 10⁶ solar masses.

At those masses, the halos can retain metals produced by high-energy pair-instability supernovae.

The simulations show that diluting the iron produced by a high-mass Population III pair-instability supernova into the gas of a roughly 10⁷-solar-mass halo naturally produces an iron abundance near the observed plateau.

The computer simulation therefore provides a physical connection between an event that occurred more than 13 billion years ago and chemical measurements made in tiny galaxies surrounding the Milky Way today.

Radiation changes the calculation

One of the most important aspects of MEGATRON is that radiation is not treated as an afterthought.

The simulation follows the buildup of radiation from early stars and its interaction with the surrounding gas.

The resulting radiation background can influence where Population III stars form and how subsequent star formation proceeds.

The four simulations consequently produce different reionization histories.

The researchers find that these differences are driven substantially by feedback-modulated escape fractions, the fraction of ionizing radiation that escapes dense star-forming environments and reaches the surrounding intergalactic medium.

This is another reason why a simple galaxy-formation model would not be sufficient for the scientific question.

The radiation field affects the chemistry.

The chemistry affects the gas.

The gas affects star formation.

Star formation determines the radiation field.

The simulation must therefore follow the feedback loop rather than calculate each component independently.

From the early universe to today’s galaxies

There is another computational trick behind the connection to present-day observations.

The researchers also evolve a matching dark-matter-only realization to z = 0.

They then use halo matching and particle-tagging techniques to follow the dynamical fate of simulated dwarf galaxies.

This allows the researchers to ask whether the tiny galaxies produced in the high redshift simulation could plausibly survive as recognizable structures in the present-day Milky Way environment.

The approach is important because the simulations themselves stop at approximately z ≈ 8.

The researchers therefore distinguish between what is directly simulated and what is inferred from subsequent dynamical modeling.

That distinction is important scientifically, and computationally.

The result is not simply a 13-billion-year movie of one galaxy.

It is a combination of high-resolution early-universe radiation-hydrodynamics and additional modeling used to connect those early structures to their potential descendants.

The supercomputer becomes a time machine

The most compelling aspect of MEGATRON may therefore be less about producing another simulated galaxy and more about creating a computational bridge between epochs that cannot be observed directly.

Astronomers cannot watch a Population III star explode inside a primordial dwarf galaxy.

Those events happened billions of years ago.

Instead, researchers can search today’s ultra-faint dwarf galaxies for the chemical fingerprints those ancient explosions may have left behind.

MEGATRON allows that hypothesis to be tested computationally.

The simulation starts with the physics of the early universe, follows the formation of the first stars, tracks radiation and chemical enrichment, generates later generations of stars, and then connects those early systems to dwarf galaxies that could survive around a Milky Way-like host.

The supercomputer effectively becomes a laboratory for experiments that nature performed only once.

Why this matters for HPC

MEGATRON illustrates a broader transition taking place in computational astrophysics.

The next generation of scientific simulations is increasingly defined not simply by larger particle counts or bigger grids, but by more complete physical models operating simultaneously at higher resolution.

Here, the challenge is the combination.

A cosmological calculation must cover a large enough region to capture the formation of a Milky Way like environment.

Within that volume, adaptive mesh refinement must resolve regions down to parsec scales.

Within those regions, the calculation must track gas dynamics, gravity, radiation, chemistry, stellar formation, and feedback.

And those calculations must be repeated under different physical assumptions to determine which predictions are robust.

That is precisely where HPC becomes a scientific instrument rather than merely a faster calculator.

The computer is enabling researchers to perform controlled experiments on physical conditions that cannot be recreated in a laboratory and cannot be directly observed in real time.

The beginning of computational archaeology

The researchers characterize MEGATRON as a vital framework for bridging high redshift simulations with the field of Galactic archaeology. This classification is highly appropriate, as the simulations suggest that specific chemical signatures observed in contemporary dwarf galaxies may function as fossilized remnants of the universe’s first stars. Furthermore, this research underscores the necessity of increasingly sophisticated supercomputing capabilities as astronomical data acquisition becomes more granular. 

Future observations from the Vera C. Rubin Observatory and next-generation extremely large telescopes will significantly expand the census of ultra-faint dwarf galaxies, providing critical datasets for validating models like MEGATRON. To maintain this scientific trajectory, computational models must evolve in tandem with observational precision. While the researchers acknowledge limitations, such as uncertainties regarding primordial star formation, nucleosynthetic yields, and the representativeness of the simulated environment, these factors do not diminish the value of the approach; rather, they define it. Each subsequent simulation serves as a refined computational experiment, narrowing the scope of plausible physical histories that could have produced the universe as we observe it today.

MEGATRON demonstrates what happens when supercomputing is pushed deep into the physics of the early universe: billions of years of cosmic history can be transformed into a numerical experiment, and the chemical fingerprints of the first stars can be tested against galaxies that still exist today.

The first stars are long gone.

Their computational fingerprints are just beginning to come into focus.

For decades, structural biology has provided scientists with high-resolution snapshots of enzymes, characterizing molecular structures frozen in specific conformations via techniques such as X-ray crystallography. However, enzymes are dynamic systems that undergo continuous conformational changes, including the opening and closing of binding pockets, side-chain rotations, and the rearrangement of water molecules and substrates. Static crystal structures inherently fail to capture these transient states; therefore, extensive computational time is required to elucidate such motions.

A recent study published in ACS Omega illustrates the significant scale of molecular-dynamics (MD) computation necessary to transform static structural data into a comprehensive portrait of enzyme behavior. Researchers investigating pyrimidine-nucleoside phosphorylase (PyNP) from Bacillus subtilis conducted an experiment involving 156 independent production trajectories, totaling 15.6 microseconds of MD simulation. 

Leveraging computational resources from the Research Center for Computational Science (RCCS) in Okazaki, Japan, alongside cloud GPU infrastructure from vast.ai, the team utilized GROMACS to analyze 13 ligands across four distinct structural states of the enzyme. Each ligand structure combination was subjected to three independent 100-nanosecond replicas, providing a robust statistical ensemble. This work underscores the transition of high-performance computing (HPC) from a mere accelerator of calculations to a critical tool for statistical sampling, ultimately allowing researchers to address a fundamental question in the field: how do enzymes function when allowed to exhibit their inherent molecular mobility?

From Four Structures to 156 Simulations

The computational campaign began with four representations of PyNP.

Three came from experimental structures, while a fourth was generated by energy minimization. They represented different positions along the enzyme’s conformational range:

• 1BRW — closed
• 5EP8min — closed-like
• 5EP8orig — semi-open
• 5OLN — open

The distinction is important because PyNP contains a flexible gate region involving residues 153–170. At its center is Tyr165, a residue positioned to move over the substrate-binding pocket as the enzyme changes between open and closed configurations.

The researchers then introduced 13 ligands into these four structural states.

That creates 52 ligand–structure combinations.

Each combination was simulated three times using independently randomized initial velocities.

The arithmetic is simple:

13 ligands × 4 conformations × 3 replicas = 156 production trajectories.

Every trajectory was 100 nanoseconds long.

Together, that produced:

156 × 100 ns = 15.6 microseconds of molecular-dynamics simulation.

Trajectories were saved every 50 picoseconds, producing 2,001 frames for each 100-nanosecond trajectory.

This is precisely the type of workload for which HPC infrastructure becomes valuable.

A single molecular-dynamics trajectory can show what happens to one molecular system under one set of initial conditions. A large ensemble allows researchers to ask whether an observed behavior persists across different molecules, conformations, and independent simulations.

The computational experiment therefore was not simply:

Run a simulation.

It was:

Run enough simulations to determine which behaviors survive statistical variation.

The Supercomputer Behind the Molecular Experiment

The production calculations used GROMACS 2025.2, a molecular-dynamics package designed for highly parallel computation.

The researchers used the AMBER ff99SB-ILDN force field for the protein, TIP3P water, and GAFF2 parameters for the ligands. The systems were solvated, neutralized, and equilibrated before production calculations were performed under NPT conditions at 300 K.

The production timestep was 2 femtoseconds, with Particle Mesh Ewald electrostatics and hydrogen-bond constraints.

At that timestep, a 100-nanosecond trajectory represents approximately 50 million integration steps.

Across 156 production trajectories, that corresponds to roughly 7.8 billion molecular-dynamics integration steps for the primary simulation campaign.

That number is not itself a measure of scientific value, but it illustrates the computational scale behind the experiment.

The researchers generated the simulations primarily on the RCCS supercomputer. For the 5EP8orig structural state, one replica was generated on RCCS while two were generated using vast.ai cloud GPUs, providing both additional computational capacity and an opportunity to examine consistency across platforms.

The study therefore represents a hybrid HPC workflow: dedicated research-supercomputing resources supplemented by cloud GPU computation.

The computational output was substantial enough that the authors deposited all 156 raw trajectories in Zenodo, along with analysis scripts.

Why 156 Trajectories Matter

The key HPC insight is that molecular simulation has a sampling problem.

An enzyme’s behavior cannot necessarily be inferred from a single trajectory. Molecular dynamics is deterministic once its initial conditions are defined, but different starting velocities can produce different microscopic histories.

That is why the study used three independent replicas for every ligand structure combination.

The researchers then analyzed the trajectories using several metrics, including:

• ligand-to-active-site distances;
• the fraction of frames in which ligands remained associated with the active site;
• residue-by-residue contact frequencies;
• MM-PB(GB)SA binding-energy estimates; and
• classifications describing ribose versus 2′-deoxyribose preferences.

This is where the HPC workload becomes scientifically meaningful.

The computer is not merely producing molecular movies.

It is producing a large statistical ensemble from which the researchers can extract patterns.

The Active Site Is Not Static

The computational results showed that the active-site pocket changes substantially between conformational states.

The calculated pocket volumes ranged from approximately 1,424 ų for 5EP8min to 2,549 ų for the open 5OLN structure.

The fully open state therefore has a pocket approximately 1.4 times larger than the closed 1BRW reference.

The authors interpret the structural progression as a transition from a contracted, substrate-trapping configuration toward an expanded, substrate-accessible configuration.

That observation establishes the computational problem.

If the pocket itself is changing size and shape, then ligand behavior cannot be understood simply by examining where a molecule sits in a single crystal structure.

The molecular machine has to be watched while it moves.

Tyr165 Emerges as a Molecular Gate

The most striking computational signal involved Tyr165.

Across the trajectory ensemble, Tyr165 showed a contact fraction of approximately 42.5% in the closed group, compared with 21.2% in the open group.

That represents a difference of 21.3 percentage points, with the reported statistical comparison producing a p-value of 1.2 × 10⁻⁵.

The physical interpretation is intuitive.

In the closed configuration, Tyr165 can sit over the active-site region, behaving like a molecular lid.

As the enzyme opens, that residue moves away from the ligand-binding region and toward solvent exposure.

The supercomputing campaign therefore converts a structural hypothesis into a trajectory-level observation:

The enzyme’s gate residue is dynamically coupled to its conformational state.

But the paper makes an important qualification.

The strongest Tyr165 signal comes primarily from 1BRW, which is a closed-state structure from the related species Geobacillus stearothermophilus, rather than from a closed-state crystal structure of the B. subtilis enzyme.

When 1BRW is excluded, the effect becomes substantially weaker.

The authors therefore do not present the result as definitive proof that Tyr165 universally behaves as a closed-state lid in B. subtilis. Instead, they identify it as a computationally supported mechanism that needs experimental confirmation using a closed-state structure from the same species.

That restraint is important.

The HPC calculation reveals a compelling molecular pattern, but the quality of the conclusion depends on the quality and comparability of the structures being sampled.

When Supercomputing Becomes a Hypothesis Generator

One of the most interesting aspects of the work is what happened after the initial 156 simulations.

The researchers used additional simulations to test whether the computationally identified mechanism could be challenged.

They performed three classes of additional all-atom molecular dynamics:

Tyr165 and other alanine mutants.

The researchers removed selected residues computationally to test predicted effects on ligand retention.

For the Y165A mutant, replacing Tyr165 with alanine substantially reduced ligand retention. The reported median late-window binding fraction was 0.02, compared with means of 0.50 for the Q153A control mutant and 0.37 for the K81A/K108A/K188A triple mutant. The reported one-sided Mann–Whitney comparison gave p = 0.049.

That is significant not because a computer has proven the biological mechanism, but because the simulation produced a falsifiable prediction.

The researchers explicitly characterize these simulations as computational support rather than experimental validation.

Adding Phosphate Changes the Question

The study also demonstrates an important limitation of classical molecular dynamics.

The initial simulations focused on enzyme–ligand complexes without the cosubstrate phosphate.

The researchers subsequently introduced an HPO₄²⁻ ion near the phosphate-binding region formed by Lys81, Lys108 and Lys188.

The phosphate interacted strongly with that lysine cluster and could approach the substrate’s anomeric carbon.

But the geometry required for the actual chemical substitution reaction was rarely observed.

The reactive in-line geometry occurred in less than 0.5% of contact frames, and in 30 of 33 systems it did not occur at all.

This result highlights a fundamental boundary between molecular dynamics and quantum chemistry.

Classical MD can model the movement and interactions of atoms using a predefined force field.

It cannot directly describe the breaking and making of chemical bonds at the electronic level.

The authors therefore point toward QM/MM and transition-state-level calculations as the next computational step.

For HPC researchers, that is an important distinction.

The computational problem does not end when the molecular dynamics run finishes.

Instead, one computational regime can identify the configurations that deserve to be examined by a more expensive method.

That is the essence of a hierarchical HPC workflow.

More Simulation Time Does Not Automatically Mean More Certainty

The study also extended selected central systems from 100 ns to 300 ns.

Those longer simulations produced retention and sugar-preference conclusions consistent with the original 100-nanosecond analysis.

But the authors appropriately qualify the result.

Only a subset of the systems was extended, and only about 46% of the original systems reached a plateau in the relevant analysis.

Consequently, the authors interpret the 100-nanosecond measurements primarily as measures of dynamic retention, rather than equilibrium affinity.

This is another important HPC lesson.

More compute is not automatically equivalent to complete sampling.

Molecular processes can occur on timescales much longer than those accessible to straightforward simulations.

Increasing a simulation from 100 to 300 nanoseconds may strengthen confidence in some observations without demonstrating that the system has reached thermodynamic equilibrium.

For the researchers, that distinction determines what the computation can legitimately claim.

From Molecular Movies to Mechanistic Maps

The study ultimately produces a more complicated picture than a simple “open versus closed” enzyme.

Different ligands interact with different combinations of active-site residues.

Gln153 emerges as an important anchor in some ligand pairs, while Lys108, Lys188, Gln153, and His82 form a broader interaction network in another case.

The computational data suggest that ribose-versus-2′-deoxyribose preferences can emerge from different residue combinations depending on the ligand and conformational state.

In other words, there is no single universal molecular switch controlling every ligand.

There is a network.

And uncovering that network is precisely where large-scale molecular simulation becomes useful.

A crystal structure can tell researchers where atoms are.

An ensemble of trajectories can begin to show how those atoms cooperate.

The HPC Pipeline Is the Experiment

Perhaps the most important lesson from the study for the supercomputing community is that the computer is not simply a supporting instrument.

The computational workflow is part of the experiment itself.

The researchers began with four molecular conformations.

They introduced 13 ligands.

They generated three independent trajectories for each combination.

They produced 156 simulations totaling 15.6 microseconds.

They analyzed thousands of molecular snapshots from those trajectories.

They identified a candidate molecular gate.

They computationally removed the gate residue.

They introduced phosphate.

They extended selected simulations.

And each stage generated another question for the next computational stage.

That is an HPC-driven scientific workflow.

The supercomputer effectively becomes a laboratory in which researchers can repeatedly perturb a molecular system, observe its response, and identify mechanisms that can subsequently be tested experimentally.

The Next Generation of the Calculation

The authors themselves identify where the computational campaign should go next.

Classical MD and endpoint MM-PB(GB)SA can capture conformational sampling and bulk electrostatic effects, but they do not resolve charge transfer, electronic polarization, or chemical bond rearrangement.

The experimentally relevant ribose-donor selectivity involves the Michaelis complex and transition state with phosphate present.

That pushes the problem toward QM/MM or QM-cluster calculations.

At the same time, understanding the kinetics of the conformational transitions will require enhanced-sampling approaches such as:

• metadynamics;
• accelerated molecular dynamics;
• transition-path sampling; and potentially
• other rare-event sampling techniques.

These approaches are often significantly more computationally intensive than standard trajectory generation, which highlights the critical role of HPC. As scientific inquiries shift from identifying visited enzyme configurations to mapping the transition pathways between them and detailing the chemical reaction mechanisms, the computational demands grow increasingly sophisticated.

The Supercomputing Takeaway

This study does not posit that high-performance computing has fully elucidated the complete mechanism of pyrimidine-nucleoside phosphorylase. Instead, it demonstrates a more significant advancement for computational science: the capacity for HPC to transform static molecular structures into statistically rigorous, sampled computational experiments.

By executing a 156-trajectory campaign, the researchers examined 13 molecular probes across four conformational states with triple replication, yielding 15.6 microseconds of molecular-dynamics data. From this ensemble, the team derived a candidate molecular lid, identified residue-level interaction networks, and generated experimentally testable hypotheses.

Concurrently, the research delineates the inherent boundaries of classical simulation: dynamic retention is not equivalent to equilibrium affinity; computational mutations do not replace empirical laboratory validation; and force-field trajectories do not capture quantum-mechanical reaction pathways. Furthermore, a statistically significant signal does not entirely negate the limitations imposed by the initial structural data.

These distinctions are not deficiencies of HPC; rather, they characterize the iterative nature of modern scientific workflows. In this model, the computer provides molecular evidence, which is then challenged by the scientist and refined through subsequent calculations or alternative computational methods. The broader potential of supercomputing in molecular science lies not merely in increased speed, but in the ability to investigate the complex mobility of molecular machines, a task that remains impossible through the analysis of static structures alone. This study, originating from four protein configurations, concludes with a compelling hypothesis regarding a single tyrosine residue and provides a clear roadmap for future, more sophisticated computational inquiry.