Quantifying selective metabolite transport for the bacterial microcompartment from Haliangium ochraceum with molecular dynamics simulations

Abstract

Bacterial microcompartments (BMCs) are protein-based organelles found in diverse bacteria that encapsulate a sequence of enzymes to accelerate specific metabolic pathways and limit toxicity by confining intermediate products. To facilitate catalysis in BMCs, reactants and products must transit across the protein shells, which are made up of multimeric protein tiles. Quantifying the permeability for transport of small molecules across the shell is a key engineering consideration to design novel catalytic microcompartments. We examine the permeability of reactants and products for two reaction pathways, the degradation of 2-aminophenol by AmnAB and the reduction of nitrite to ammonium by NrfA, through a computational lens to quantify permeability for these metabolites at the nanoscale. From Hamiltonian replica exchange umbrella sampling simulations, we determine that the energetic barriers to permeation for the relevant metabolites along both pathways are small, with permeability coefficients in the range of 0.1–4 cm s⁻¹. The high permeabilities calculated for each metabolite through the inhomogeneous solubility diffusion model are corroborated by enzyme activity measurements in vitro that point to enzymatic activity within the shell system. We expand on these findings by quantifying the scale of transport and catalysis that can be supported by in vitro BMC systems. The results highlight the importance of substrate permeability across BMCs within a biological context, and represent steps toward generating synthetic shells or bioengineering novel nanoreactors.

---

Introduction

Bacterial microcompartments (BMCs) are found across multiple bacterial phyla,¹ encapsulating enzymes and thus metabolic reactions within a protein shell.^2,3 Enzyme encapsulation within BMCs is thought to enhance survival by facilitating reaction coupling along a pathway, creating locally high concentrations of reaction intermediates.^4,5 The high intermediate product concentration can then be consumed by other co-confined enzymes in the reaction pathway, with both catabolic and anabolic reaction pathways found in nature facilitating diverse chemistry.^1,4,5 Crucially, compartmentalization is suggested to create a barrier to diffusion for products and reactants that enhance catalytic activity.^3,6,7

The diffusive barrier formed by the shell must be large enough to retain a high intermediate product concentration, while simultaneously only forming a small barrier to permeation for the ultimate inputs and outputs for the pathway. Since BMC shells are made from multiple proteins that feature a central pore,⁸ one possibility is that these pores in the center of hexameric, pentameric, or trimeric subunits regulate metabolite flux by acting as a semi-permeable barrier, increasing the permeability of products and reactants, while lowering the permeability of intermediate products that may be harmful to the cell.^6,9 Different electrostatic environments across the shell may facilitate this sorting,¹⁰ and prior molecular simulations have identified barriers to transit across BMC pores depending on permeant characteristics like charge and size.^11–16 Past studies focused on pores coupled to their native substrates. However, for BMC architectures to live up to their engineering potential as nanofoundries for specialized products,^3,17 permeation for non-native substrates also must be considered, particularly where evolutionary pressures have not optimized permeabilities for specific reactants or products.

Through this engineering lens, we interrogate the partial encapsulation of two distinct metabolic pathways, the reduction of nitrite to ammonium by the cytochrome c nitrite reductase (NrfA) enzyme¹⁸ and the catabolism of 2-aminophenol (2-AP) into picolinic acid and further downstream products¹⁹ (Fig. 1). NrfA catalyzes a fundamental reaction in the global nitrogen cycle^18,20 and is studied extensively in part as a route to sustainable nitrogen cycling for agriculture. Bacterial degradation of aromatic substrates such as 2-AP towards central metabolism is essential for growth on mixed waste streams, and bacteria with these enzymes are used to for their remediation.²¹ 2-AP 1,6-dioxygenase (AmnAB), the enzyme that commits 2-AP towards degradation, has an aldehyde intermediate that spontaneously cyclizes. Since the intermediate aldehyde can be reactive, the pathway is confined within BMCs in Acintobacteria,¹⁹ and so we are translating this pathway into a different BMC context.

Fig. 1 Encapsulated pathways under study. (A) NrfA converts nitrite into ammonium when provided reducing equivalents, which in this study is dithionite added exogenously. (B) AmnAB is the first step of a larger pathway in 2-aminophenol (2-AP) degradation, creating an intermediate aldehyde that cyclizes into picolinate. (C) Presents a schematic BMC illustrating the fundamental permeation problem created by encapsulating specific enzymes.

To feed catalytic pathways within BMCs, reactants and products must permeate across the shell. All BMCs possess a proteinaceous shell composed of multiple subunits (BMC-H, BMC-T, and BMC-P; Fig. 1C), each contributing distinct flux properties.²² Molecular and genetic studies indicate that reactants, products, and co-factors for the reaction enter BMCs through hexameric BMC-H and trimeric BMC-T shell proteins that separate the BMC interior from the bacterial cytoplasm.^23,24 Molecular and genetic studies indicate that reactants, products, and cofactors for the reaction enter BMCs through pores in hexameric (BMC-H) and trimeric (BMC-T) shell proteins that separate the BMC interior from the bacterial cytoplasm.^23,24 The key questions here are whether the NrfA and AmnAB enzymes can be encapsulated in non-native systems, such as the BMC from Haliangium ochraceum with well-established molecular toolsets for engineering and BMC design,^25,26 and whether the permeability of these BMCs is sufficient for robust enzyme catalysis within the confined BMC shell.

Quantifying permeability experimentally is difficult. Many particle tracking tools would substantially change metabolite dynamics, and radiolabling approaches may not be compatible with the timescales needed to track BMC permeation. The simultaneous spatial and temporal resolution of molecular simulation obviates many of these difficulties, as individual particles can be tracked as they transit across BMC components.¹¹ When combined with enhanced sampling techniques and permeability theory,²⁷ molecular simulation is uniquely suited to directly quantify permeability for the small molecules listed in Fig. 1 across both BMC-H and BMC-T shell components. Together with complementary experiments for assembled shells, we find that the molecular permeability is not rate limiting for these particular reaction substrates when confined within the H. ochraceum BMC (HO-BMC).

Methods

This study has two complementary parts. In the first, we determined the permeability coefficients for the small molecule metabolites introduced in Fig. 1 through molecular simulation to establish an estimate for the permeability. We then corroborate the high permeability coefficients obtained via simulations through in vitro work to track activity for these pathways once encapsulated.

Determining metabolite permeability through molecular simulation

As we do not know a priori how fast metabolites may permeate the shell components, our approach was to take the flattened sheet structure from prior simulations of the HO-BMC²⁸ as the basis for simulation. This starting point was used to prepare Hamiltonian replica exchange umbrella sampling (REUS)²⁹ to overcome potential limitations due to limited accessible simulation timescales. This approach has been used to determine passive membrane permeability through the inhomogeneous solubility–diffusion model (ISDM), which was developed to measure permeability in heterogeneous lipid bilayers and is now considered the gold standard for calculating membrane permeability coefficients.^27,30 The ISDM is derived from the steady-state flux and assumes equilibrium across the membrane, allowing spatial variation of both the diffusion coefficient and the free energy profile within the membrane interior. The potential of mean force (PMF), denoted as ΔG(ξ), and the position-dependent diffusion coefficient, D(ξ), can be obtained from MD simulation data and are mathematically related to the resistivity (R) and permeability (P) as described in eqn (1).³¹

In eqn (1), β refers to the beta thermodynamic term, defined as (β = 1/k_BT), and ξ represents the collective variable describing the relative position of the metabolite along the BMC pore axis. The integration limits, ξ_l and ξ_u, correspond to the lower and upper bounds of the pore axis, respectively.

While clearly different from a lipid bilayer, the BMC shell is also a “membrane”, and so the same ISDM formalism can be used to determine permeability. With the ultimate goal of determining permeability coefficients in mind, we prepared eight independent simulation systems, where a small planar BMC shell fragment containing three hexamer and one trimer tile has had five copies of a single metabolite placed along the axis normal to the pores present in the BMC. In this arrangement, one molecule was placed along each hexamer pore, and two molecules were used to sample along the double-stacked T2 or T3 trimer unit used in prior HO-BMC models²⁸ (Fig. 2A). Thus, each of these metabolites is placed along the z-axis to sample the passage of the individual metabolites across each shell protein (BMC-H or BMC-T). Since we needed to sample along this axis fully, we conducted steered molecular dynamic (SMD) simulations to drive the metabolites across each BMC pore to seed the REUS starting points. The details of this approach, and the analysis that follow, are given in the subsequent subsections and SI.

Fig. 2 Pictorial representation of metabolites traversing the stacked trimeric and hexameric pores. Metabolites are shown in yellow. Hexameric and trimeric units are colored magenta and purple, respectively. The arrows represent the collective variable, which spans from ±60 Å for stacked trimers and ±30 Å for hexamers. Panels (A) and (B) denote the side and top view for a metabolite crossing through the pore, with panel (B) emphasizing that the metabolite was confined to sample within 15 Å of the pore center. Panel (C) illustrates the BMC-hexamer and stacked BMC-trimer along with the metabolites shown in the solvated simulation box to get a sense of the scale. Panel (D) illustrates the periodic boundary representation of the asymmetric unit, with the 3 hexamers and 1 trimer colored while the periodic images are ghostly gray.

System preparation and biased simulations

The system construction follows the approach from Raza et al.²⁸ using a small planar shell fragment as the starting point. All the systems were solvated with TIP3P water model³² and the charges were neutralized by adding Na⁺ and Cl⁻ ions such that the resulting concentration is around 0.15 M. Subsequently, SMD simulations were performed in NAMD 3.0a9³³ or newer (Fig. S1). Because we were interested in the metabolite transit across the BMC pore, the natural choice of the pulling variable was the position along the Z axis that was co-linear with the pore (Fig. S1). The spring constant used for the initial pulling simulations was determined by performing different trial simulations across a range of spring constants, leading to 30 ns SMD pulls with a 5 kcal mol⁻¹ Å⁻² force constant.

For pulling, five metabolites were placed, three in line with the hexamer pores and placed 30 Å from the membrane center and two metabolites were placed on each side of the trimer dimer at 60 Å and −60 Å (Fig. 2A and B). The SMD simulations were performed in two steps. In the first step, for hexamers the metabolites were pulled from 30 to −30 Å, whereas one trimer side metabolite was pulled from solvent to the center (60 to 0 Å), while the other was held at −60 Å (Fig. S1). In the second step the hexameric side metabolite were pulled from −30 to 30 Å, while the molecule pair pulled through the trimer were moved from −60 to 0 Å and 0 to 60 Å, respectively (Fig. S1). The collective variables module in NAMD³⁴ was further used to define the cylindrical coordinate system specific to the metabolites of interest so that each metabolite was confined to a specific pore. A flat-bottomed potential with a force constant of 5 kcal (mol Ų)⁻¹ restrained each small molecules from moving more than 15 Å away from pore axis.

For REUS simulations, 128 independent replicates were run simultaneously, splitting the 60 Å distance each metabolite needs to sample along the reaction coordinate into approximately 0.5 Å wide windows. The initial seeds for the 128 replicas were taken from equally spaced timesteps during the second step of the SMD. As was done during the SMD, the force constant for the harmonic restraint that confines the metabolite within the window along the pore axis was 5 kcal mol⁻¹ Å⁻².

In prior simulations, the hexameric or trimeric units would tilt slightly during simulation,²⁸ which would have negative consequences on the z-coordinate based reaction coordinate being sampled for REUS. To eliminate tilt, we applied an additional restraint to prevent tilting for the individual hexamer or trimer tiles within the surface. The tilt restraint was calculated on a tile by tile basis using the center of mass of two groups of atoms to define an axis that should be aligned with the z-axis. The force constant applied was very high, with a k = 10 000 000 kcal mol⁻¹ scaling factor applied to the cosine of the angle between the defined axis and the z-axis. In practice, a 1 degree tilt angle under these circumstances would contribute 0.2 kcal mol⁻¹ to the overall energy (k(1 − cos(1°))²), while larger tilts would be substantially penalized.

Exchanges between adjacent replicas during REUS were carried out every picosecond, with exchange ratios near the 20% optimum.³⁵ The final duration for each REUS simulation is listed in Table S1, with metabolites demonstrating converged free energy profiles at different levels of sampling (Fig. S2 and S3). Most of the metabolities were converged within 80 ns of simulation time. Dithionite was difficult to converge and required extended simulations up to 150 ns.

For both the REUS and SMD simulations, we used versions of NAMD, either 2.14 for minimization, or NAMD version 3.0b5.³³ The CHARMM36m force field was used for the protein component.³⁶ The systems were solvated using the TIP3 water³² model. The metabolites were parameterized initially through CGenFF,³⁷ with additional optimization for problematic groups using the FFParam tool kit.³⁸ The details of parameterization are given in the SI (Fig. S4, S5 and Tables S2–S5).

A Langevin barostat and thermostat were used to maintain the pressure and temperature at 1 atm and 298 K, respectively, with 1 ps⁻¹ damping.³⁹ Hydrogen bonds were handled with the SETTLE algorithm to enable 2 fs timesteps.⁴⁰ The long-range non-bonded Lennard Jones (LJ) cutoff was 12 Å. Long-range electrostatic interactions were calculated with the particle mesh Ewald (PME) grid with 1.2 Å.^41,42 The switching between electrostatic and non-bonded interactions were done after 10 Å. LJ correction was applied to conserve energy during switching.⁴³ Energy minimization of the system was initially performed using the 1000 steps of conjugate gradient in NAMD. The system was then allowed to equilibrate for 50 ps in NPT ensemble using a 5 Å margin to allow the box to adjust after any distortion from minimization prior to transitioning to the GPU-resident integrator during the SMD step.

Simulation analysis and permeability determination

Analysis of the simulation trajectories was carried out using the combination of built-in and purpose built VMD⁴⁴ python scripts and features, leveraging the numpy and scipy libraries^45,46 to handle numerically complicated tasks.

To determine a molecular-scale permeability via eqn (1), we needed both a free energy and diffusivity profile along the relevant reaction coordinate, in this case the z-axis co-linear with the pores within the BMC-H and stacked BMC-T subunits. The estimation of the free energy profile was performed via BayesWHAM⁴⁷ directly from the collective variable trajectories generated during the REUS simulations. The local diffusivity was estimated from the variance and time autocorrelation along the reaction coordinate.^48–50 An exponential fit to the autocorrelation function at time scales smaller than the exchange frequency was used to calculate the temporal autocorrelation in order to account for the exchanges introduced by REUS sampling. Using the approach from eqn (1), the computed profiles for diffusivity and free energy are then fit to a smooth spline and integrated using Simpson's rule to find the permeability.

We deliberately computed the permeability across the pore, as we restrained the metabolite during REUS to within a cylinder with a 30 Å diameter (Fig. S6). This was essential to maintain the metabolite position, preventing it from transiting the shell outside of the pore, thus improving sampling near the BMC surface. However, the BMC-T and BMC-H tiles are physically larger than the cylinder the metabolites were constrained within. Assuming that the permeability is zero outside of the pore we sampled, and that we accurately modeled the permeability within this confined space, we estimated an effective permeability (P_effective). P_effective is the permeability along the hexameric (P_H) and trimeric (P_T) pores, multiplied by the fraction of confined cylinder area with radius (r = 15 Å) (eqn (2)).

In eqn (2), the A is the total area of the simulation cell in x–y plane. This area covered all four pores (3 BMC-H and 1 BMC-T pore) present in the shell. Since all the pores were of roughly same size, the key scaling factor is . In all the simulations, the confined radius was 15 Å, and the surface area of the box was ≈16 200 Ų, leading to the factor of 0.0437. Upon incorporating the contribution of all four pores in the shell, the effective permeability across the total surface area was reduced to approximately 16% of the permeability estimated using the cylindrical restraint configuration.

Quantitative analysis of protein–metabolite interactions

To quantify protein–metabolite interactions, we evaluated the number of metabolite approaches to the protein within a specified distance cutoff. As contact numbers are highly sensitive to choice of distance cutoff, we adopted an exponentially weighted contact function that enables the inclusion of fractional contributions from interactions occurring near the cutoff boundary. This approach provided a smoother and more physically meaningful representation of contact behavior, and has been employed in prior studies to mitigate the limitations of binary contact definitions.^51–53

In eqn (3), the contact measure only depends on d_ij, the distance between atoms i and j.

Experimental methods

Beyond the simulations, in vitro experiments were used to demonstrate confined catalysis for the AmnAB and NrfA within a BMC from H. ochraceum. The NrfA and AmnAB experiments were designed as straightforward proofs of concept to demonstrate the plausibility of the calculated permeabilities. As such, these experiments were intentionally kept simple, are not meant to define the amount of encapsulated enzyme or provide quantitative rates, and also do not yet reflect fully optimized reaction conditions.

NrfA loaded BMC shells were prepared using a chaotropic in vitro assembly protocol including all 6× His tagged shell components.⁵⁴ Shell components were mixed in ratios in a solution of BMC-H tiles and urea to form full shells. NrfA from Geobacter lovleyi was inserted onto a pBAD202/D-TOPO vector (Invitrogen) with the addition of a SpyC001 on the N-terminus (see details in SI, Fig. S7). SpyC-NrfA was then induced with 0.02% arabinose and expressed in Shewanella oneidensis at 30 °C for 16 hours.^54,55

SpyT-BMC-T1-6× His that had incubated for 1 h with G. lovleyi SpyC-NrfA-Strep-tag II was added last to initiate BMC assembly. Complete shells containing NrfA were purified from the assembly mixture by passage through a Cytiva S200 size exclusion chromatography column.

To further ensure minimal exposed Strep-tagged Spy-catcher-NrfA in these samples, shells containing NrfA were passed through a gravity column with IBA Strep-Tactin XT 4Flow resin. The flowthrough was collected and concentrated on a Millipore 3 kDa molecular weight cutoff Amicon Ultra centrifugal filter. Before assaying, extra BMC-P was added to samples to a concentration of 0.47 mg mL⁻¹ of monomer to further ensure capping of the BMC vertices. Encapsulated NrfA activity was assayed using reduced methyl viologen where the decrease of 600 nm absorbance was used to measure the oxidation of methyl viologen as NrfA reduced nitrite. Samples were assayed in technical triplicate on a BioTek Synergy HTX plate reader in a Coy anaerobic chamber.

To obtain AmnAB loaded BMC shells, we used an in vivo expression approach with two compatible plasmids. The first plasmid, which encodes AmnA¹⁹ fused to the BMC-T1 protein from H. ochraceum⁵⁶ and also carries genes for BMC-T2, BMC-T3, BMC-H, and BMC-P, was co-harbored with a second plasmid expressing AmnB (see details in SI, Fig. S8). Protein expression in E. coli was induced with 0.1 mM isopropyl β-D-1-thiogalactopyranoside.

The cells were disrupted by using a French pressure cell. Cargo-loaded BMCs were purified by using a 30% sucrose cushion and Mono Q column chromatography according to previously described procedures.²⁶ Non-encapsulated AmnAB also was purified by methods described previously.¹⁹ The 2-AP 1,6-dioxygenase activity for 2-AP was monitored by using a UV-visible spectroscopic assay at 380 nm,¹⁹ the absorbance maximum of the intermediate 2-aminomuconic acid-6-semialdehyde, whereas there is no significant absorbance at this wavelength for the product and reactant. A total of 200 μL of purified AmnAB-loaded BMCs (0.3 mg mL⁻¹) was incubated with 50 μL of 2-aminophenol at various concentrations to measure enzymatic activity.

Results

Free energy and local diffusivity

The key goal for this work was to determine the permeability for metabolites involved in NrfA or AmnAB catalysis, and to do so we used the ISDM shown in eqn (1). This analysis required profiles along our chosen reaction coordinate for the free energy and diffusion of small molecules across hexameric and trimeric pores based on MD simulations using REUS to integrate this equation. The profiles shown in Fig. 3 highlight the range of energetic barriers possible with different metabolites.

Fig. 3 Metabolite free energies across the stacked trimeric (A) and hexameric (B) pores. The free energies for the 2-AP pathway metabolites are shown in the top panels, with a blue background matching Fig. 1, while the nitrite to ammonia pathway is in the bottom panel with a green background. The reaction coordinate for both hexamers and trimers considers the center of mass of the hexamer or trimer pore as the z = 0 position, and the free energy in solution is at the left and right hand edges. These profiles neglect the first 20 ns of simulation data, which we consider as the equilibration period after pulling in the SMD. The solid lines indicate the free energy profiles calculated by BayesWHAM,⁴⁷ while the dashed lines are the spline-fit integrated by eqn (1).

Free energy profiles for permeation across BMC tiles

Typically, the free energy profiles were negative within the BMC shell, indicating that the small molecules interacted with the protein rather than remaining in solution. The largest positive free energy for any small molecule involved with the 2-AP pathway was below 1 kcal mol⁻¹, while small molecules within the BMC experienced free energy minima between −2 kcal mol⁻¹ and −3 kcal mol⁻¹. The negatively charged picolinate was perhaps the most unique small molecule within this set, with the most favorable interaction with the interior of the trimer, up to approximately −3 kcal mol⁻¹. The uncertainties for small molecule free energies within the trimer pore were larger than those for the equivalent hexamer. Based on our simulation design, the hexamer profiles combined sampling through three independent pores to generate the free energy estimates, while the trimer estimates depend on sampling through the single pore (Fig. S9). Regardless, the relatively small peaks and valleys for AmnAB metabolites show how little a barrier the shell poses to these metabolites.

The profiles associated with the NrfA pathway were largely similar to those for AmnAB. The negatively charged molecules nitrite and dithionite have multiple peaks and valleys, but generally had minimal barriers to permeation through the trimer. Dithionite transport through the hexamer had a notable free energy barrier of around 2.5 kcal mol⁻¹, which was not observed in the larger trimer pore. It is unlikely that dithionite is close to the limit for the size at which small molecules can permeate across the shell. Instead, the two negative charges on opposite ends of the molecule may prevent the small molecule from accessing specific regions of the reaction coordinate by mediating strong interactions along the pore. On the other hand, the positively charged ammonium had the largest free energy barrier of any molecule tested, with barriers near 4 kcal mol⁻¹. This in itself was not surprising, as positively charged ions also had large barriers when studied in related carboxysome systems.¹²

Diffusivity profiles for permeation across BMC tiles

Together with the free energy profile, local diffusion coefficients for all the metabolites are essential to evaluating permeability, and are shown in Fig. 4. All cases exhibited the same qualitative behavior; diffusion was largest in the aqueous region, then it dropped as permeants approached the pore, with a minimum value typically close to where bottlenecks are located. The stacked trimer contains three internal aqueous pockets-one per trimer and an additional interfacial aqueous region (Fig. S10). The minimum diffusion coefficient was at the two gates between these interior compartments. The hexamers tended to have a single, poorly-defined minimum, as the bottleneck there was simply the constriction where the six monomers within the hexamer make their closest approach.

Fig. 4 Diffusion of metabolites across the trimeric (A) and hexameric (B) pores. The solid lines indicate the position-dependent diffusion values, whereas the dashed lines represent the spline-interpolated fits that were integrated to compute permeability (eqn (1)). The first 20 ns of the data was considered as equilibration period.

The most interesting individual values for diffusion were seen in the NrfA pathway metabolites. Quantitatively, for both the hexameric and trimeric pores, dithionite exhibited substantially slower diffusion. We hypothesize that despite its small size, dithionite's relatively high molecular weight due to its sulfur atoms being much heavier than carbon or oxygen contributes to this reduced diffusivity. Nitrite and ammonium ions exhibited similar diffusion profiles across the pores, despite the disparate charges for these molecules and varying molecular weights (Fig. 2).

In the case of metabolites related to 2-AP, there were no significant differences between the 2-aminomuconic acid-6-semialdehyde and 2-aminomuconate-6-semialdehyde molecules while passing through the hexameric or trimeric pores. These two molecules differ only in their charge, and are neutral or negatively charged, respectively, which would suggest equal diffusion regardless of charge. However, comparing between picolinic acid and picolinate, the diffusion of the picolinate was slower compared to picolinic acid, showing altered behavior compared with other biomolecules. 2-AP is quantitatively similar to other 2-AP pathway molecules in both pores.

Overall, while the diffusion coefficients were inherently noisy, the smoothed profiles were quite similar. Across the full permeation pathway, each metabolite exhibited relatively modest variability in diffusivity; for any given molecule, the maximum difference between the highest and lowest diffusivity values was less than a twofold change. Due to the structure of eqn (1), where diffusivity is directly related to permeability while the free energy profile is exponentially connected to the permeability, these small changes in the diffusivity were relatively minor contributors to the overall free energy profile. Instead, what is most important to realize here is that the diffusivity for these small molecules was on the order of around 80 Ų ns⁻¹, which sets the relative rate at which these molecules permeate across the shell. This basal rate was modified by the exponential factor from the free energy profiles shown in Fig. 3.

Metabolite permeability through the BMC pores

Together, the free energy profiles and the local diffusion constants can be combined via the ISDM (eqn (1)) to estimate the permeability coefficients (Table 1). Overall, the range for the permeability coefficients was rather small, and are similar in magnitude to estimates for other small molecules using the same methodology.¹⁶ All permeability coefficients (P_effective) reported in Table 1 are separated by less than an order of magnitude, yielding permeabilities that can be easily graphed on a single axis (Fig. 5). For comparison, for diverse aromatic metabolites across lipid bilayers the permeabilities can range by 6 orders of magnitude.⁵⁷ This suggests that molecular permeabilities for molecules as diverse as what are shown in Fig. 1 are all quite similar as they go across a BMC.

Fig. 5 Permeability of metabolites across the hexamer (P_H) and trimeric (P_T) pores. The cyan and green color background indicates the 2-AP and NrfA pathways, respectively, as shown in Fig. 1.

Table 1 Permeability (cm s⁻¹) values for metabolites across the hexameric (P_H) and trimeric (P_T) pores. The effective permeability (P_effective) was calculated to adjust the cylindrical confinement of molecules crossing the pore eqn (2). Individual table rows are colored according to Fig. 1

Compound name	PH	PT	Peffective
2-Aminophenol	23.1 ± 0.1	10.8 ± 0.1	3.5 ± 0.1
2-Aminomuconic acid-6-semialdehyde	11.1 ± 0.1	6.5 ± 0.1	1.7 ± 0.1
2-Aminomuconate-6-semialdehyde	11.0 ± 0.1	11.4 ± 0.1	1.9 ± 0.2
Picolinic acid	15.4 ± 0.1	7.3 ± 0.1	2.3 ± 0.1
Picolinate	13.2 ± 0.1	9.1 ± 0.1	2.1 ± 0.1
Ammonium	0.3 ± 0.0	0.6 ± 0.0	0.1 ± 0.0
Nitrite	14.2 ± 0.1	12.2 ± 0.2	2.4 ± 0.2
Dithionite	2.3 ± 0.0	8.6 ± 0.1	0.7 ± 0.1

---

While the overall permeabilities were quite similar, there were some trends. The two slowest molecules were ammonium, which is positively charged, and dithionite, which is unusually heavy and doubly negatively charged compared with other molecules in this set. The fastest molecule was distinct between hexamers and trimers, and were relatively close together compared with the substantial multiplicative differences between the slowest permeants. Hexamer and trimer pores, despite different geometries for their bottleneck radii,^28,56 had similar permeabilities. The smaller hexamer pore appeared to have virtually indistinguishable permeation rates when compared to larger trimer pores, with the permeability order flip flopping multiple times in Fig. 5. In the NrfA pathway the negatively charged dithionite was more permeable through the trimer as compared to the hexameric pore (Fig. 5), while nitrite was the reverse. In contrast the positively charged ammonium had almost similar permeability across both the pores. The 2-AP pathway permeabilities were similarly scattered, with the hexamer and trimer exchanging which path was faster.

The finding that the trimers and hexamers had similar permeabilities would seem to be counterintuitive with the observations in Fig. 3, where the hexamers tended to have shallower free energy profiles than the trimers. However, eqn (1) explicitly integrates over the length of the reaction coordinate. Since our trimers are a double stack, the longer path length within the trimer more or less cancels out the more favorable interactions along the trimer path due to the larger bottleneck radius than the equivalent hexamer. We would thus anticipate that single stacked T1 trimers would exhibit faster permeabilities.

The reason we have multiple columns in Table 1 was due to the geometry of a real BMC. The first two columns of Table 1 are what we would compute directly from our simulations. However, as those simulations restrained the metabolites along the pores, and there was additional surface area that was inaccessible due to the applied constraints. Since that extra surface area is impermeable to metabolites as there is no pore present, the effective permeability (Fig. S11) will be smaller than what we report based strictly on geometry (eqn (2)). This correction was important, as it refined the calculated permeability values that initially overestimated rates to an estimate of what would be observed across a whole BMC shell.

Permeability coefficients in the range from 1–10 cm s⁻¹ like we see in Table 1 are very fast compared to other biological transport processes. Lipophilic compounds through a lipid bilayer have permeability coefficients in this range.^57–59 Hydrophilic compound permeability coefficients are substantially slower across a bilayer, with permeation coefficients for water approximately 1 million times smaller.⁶⁰ Even fast channels like aquaporins have permeability coefficients that are 100 times smaller than what we see here for BMC permeabilities.^61,62 The chief difference is that natural systems do not array aquaporins adjacent to one another, so the pores are relatively far apart compared to the tightly packed BMC tiles. Thus the high pore density in BMC shells suggests that BMCs are evolutionarily optimized for high flux, as a hypothetical “closed” tile without a pore has not yet been identified in a BMC.

Electrostatic charge drives molecular transit across the pore

From the free energy and diffusivity profiles, perhaps the biggest surprise was the difference in behavior between picolinic acid and picolinate. Specifically, the free energy differences between picolinic acid and picolinate shown in Fig. 3 were quite large given the rather small change in structure, especially compared with 2-aminomuconate-6-semialdehyde or 2-aminomuconic acid-6-semialdehyde, which also only differ by a proton. We take the time here to evaluate what exactly the differences are in the free energy profile and how these relate to observed interactions within these simulations.

The free energy profile for picolinic acid in the trimer indicated two favorable regions (ΔG ∼ −3.0 kcal mol⁻¹) around the interface of stacked trimers and close to the concave side of the upper trimer (Fig. S9). In contrast, the free energy profile for picolinate was smooth and associated with multiple low energy minima (ΔG ranging from 0 to 1 kcal mol⁻¹), with the lowest point occurring at around the center of trimer. Fig. 6 shows the residue interaction network by picolinate (Fig. 6A) and picolinic acid (Fig. 6B) near this position. To quantify these interactions more closely, we turned to the contacts made between the metabolite and different residues in the protein via eqn (3). The normalized contact analysis between picolinic acid and the trimeric assembly highlights strong interactions with Arg27, Arg70, and Tyr172, as reflected by distinct peaks in Fig. S12. Picolinic acid made contact with the trimer in basically the same way, with peaks in the same locations, but much shallower, suggesting that the uncharged molecule made fewer contacts with the protein as it transits the pore. The trimer pore region is enriched in basic amino acids (Fig. S13 and S14), which interact more favorably with the negatively charged picolinate compared with picolinic acid. This electrostatic complementarity likely accounts for the pronounced peak observed in the contact profile (Fig. S12), as well as the comparatively deeper PMF.

Fig. 6 Interaction of metabolites with the shell. Interaction formed between picolinic acid (A) and picolinate (B) in the trimers. Panel (C) and (D) plots the interaction of dithionite in hexamers (C) and trimers (D), respectively. The snapshot shows amino acids present within 5 Å of the metabolites. The contacts are quantified in Fig. S12 and S15. For clarity only residues in a plane close to the metabolites are shown.

Surprisingly, for picolinate, the lowest free energy minimum was observed at approximately −5 Å. The contact analysis (Fig. S12) revealed multiple interaction peaks corresponding to several charged residues (Arg27 and Arg70) and aromatic residues (Tyr172). The free energy plot of picolinate in the trimer indicated two favorable regions (ΔG ∼ −3.0 kcal mol⁻¹) around the interface of stacked trimers and close to the concave side of the upper trimer (Fig. S9). In contrast, the free energy profile for picolinic acid is smooth and associated with a multiple low energy minima (ΔG ranging from 0 to 1 kcal mol⁻¹), with the lowest point occurring at around the center of trimer. Fig. 6 shows the residue interaction network formed by picolinate and picolinic acid near this position.

In contrast, the dithionite molecule had a lower free energy minimum as it crossed the trimer (ΔG ≃ 6 kcal mol⁻¹, Fig. 6D), indicating a smooth passage from the stacked trimers. While dithionite showed a barrier for permeation near the pore bottleneck as it crosses the hexamer (Fig. 6C), it also had two minima near −12 Å and 8 Å where dithionite would interact favorably with the hexamer. From Fig. 6C and D we see that interactions with charged amino acids such as Glu and Arg were responsible for reducing the free energy in these regions, facilitating dithionite transfer across the shell. Correspondingly, contact analysis plots (Fig. S15) exhibited pronounced peaks associated with Glu and Arg residues. Similar trends were also observed for the nitrite, 2-aminomuconic acid-6-semialdehyde, 2-aminomuconate-6-semialdehyde, and ammonium for the hexamer and trimer (Fig. S16 and S17). Collectively, these findings demonstrate that the presence of energetically favorable electrostatic interactions-primarily involving charged residues within and adjacent to the pore region-are what dictate the free energy barriers reported in Fig. 3, thereby dictating the efficiency and likelihood of metabolite translocation across the BMC shell.

Experimental validation of metabolite permeability in encapsulated systems

The simulation findings suggest that the permeability for both of the pathways encapsulated within a BMC will be high, and unlikely to be a major bottleneck to allowing catalysis to happen within a BMC. As a proof of concept, we also carried out experiments using established in vivo assembly methods and newly developed in vitro assembly techniques⁵⁴ to test whether both reactions shown in Fig. 1 are indeed feasible within a minimal BMC shell. Our expectation was that since permeation rates were very high, the reactions from Fig. 1 should proceed unimpeded, even when the critical enzymes are encapsulated.

NrfA encapsulation and methyl viologen-based activity assay

To evaluate the catalytic efficiency of encapsulated NrfA, the simplest test is an enzyme activity assay using Spy-tagged hexamer, trimer, and pentamer (HTP) shells loaded with complementary SpyC-NrfA. Only encapsulated NrfA enzymes will be present after purification. NrfA catalyzes the six-electron reduction of nitrite (NO₂⁻) to ammonium (NH₄⁺), as illustrated in Fig. 1, and so if it is encapsulated, nitrite and ammonium must both permeate to allow this reaction to proceed. Because neither nitrite nor ammonium are particularly good reporters, we used methyl viologen as an artificial electron donor, and tracked its oxidation state. Reduced methyl viologen is blue, but when it is oxidized the color visibly changes, leading to a loss of absorbance at 600 nm.

As shown in Fig. 7, the absorbance rapidly decreased when nitrite was added to the capped BMC HTP shells that contain NrfA. When NrfA was missing or no nitrite was added, the 600 nm absorbance was constant. With increasing concentrations of nitrite (either 500 μM or 1000 μM), however, the absorbance decayed more quickly as methyl viologen was oxidized. The timescale for a full reaction of approximately 2 min was also similar to prior one-step reactions performed inside of a HO shell.^7,54

Fig. 7 HTP shells with cytochrome c nitrite reductase (NrfA) cargo exhibit activity. Enzyme activity was measured via methyl viologen activity assay. The loss of absorbance of 600 nm indicate methyl viologen oxidation which is coupled to nitrite reduction by NrfA. Capped HTP shells with NrfA cargo oxidized methyl viologen. All the experiments were repeated in triplicate. The plot illustrates the mean and standard error of the recorded results.

Spectroscopic assessment of catalytic activity for encapsulated AmnAB

Evaluating the encapsulated AmnAB enzyme that results after purification was done in parallel to the NrfA work. AmnAB catalyzes the oxidative ring cleavage of 2-AP, producing 2-aminomuconic acid-6-semialdehyde, which subsequently undergoes spontaneous cyclization to form picolinic acid.¹⁹ The intermediate 2-aminomuconic acid-6-semialdehyde exhibits a characteristic absorbance peak at 380 nm (Fig. 8), which is different from either the reactant 2-aminophenol, or the ultimate picolinate/picolininic acid product. Thus, for an encapsulated system, we would expect a transient increase in absorbance at 380 nm shortly after starting the reaction, but quickly reducing in absorbance as the reaction proceeds to the product that does not absorb at 380 nm. This is exactly what we see in Fig. 8, suggesting that encapsulating AmnAB is not detrimental to its function, and that substrates and products can readily permeate just as was calculated from simulation.

Fig. 8 2-Aminophenol-1,6-dioxygenase activity of (A) encapsulated or (B) free AmnAB as monitored by the absorbance of transiently produced 2-aminomuconic acid-6-semialdehyde. BMCs containing AmnAB were constructed and purified as described in the methods section and used at 0.3 mg ml⁻¹ total protein. Free AmnAB, purified as described¹⁹ was used at 0.5 mg ml⁻¹. The samples were incubated with 10 or 100 mM 2-aminophenol, and the absorbance at 380 nm was monitored using a plate reader.

Discussion

We see from simulation that estimated permeabilities range between 0.1 to 3.5 cm s⁻¹ (Table 1). Permeabilities in this range have been reported for other BMC systems using the same approach,¹⁶ and are 100–1000x faster than permeabilities determined at high metabolite concentration through unbiased simulations.^11,16 Crucially, BMC permeabilities are at least a million times higher than permeabilities would be for lipid bilayers for the charged molecules we are interested in here, with charged metabolites often having permeabilities in the range of 10⁻⁶ cm s⁻¹.⁵⁷ At this high level of permeability through a BMC shell, small molecule permeation is unlikely to be the rate-limiting factor for catalysis in vitro. Instead, the overall catalytic efficiency is more likely governed by intrinsic enzyme kinetics and diffusion limitations within the confined shell environment. The key challenge, therefore, lies in determining how far we can drive the conversion of reactants to products before encountering these enzyme-specific kinetic and diffusion constraints inherent to encapsulated enzymes.

To test these limitations, we made a few assumptions to bound parameters for a theoretical model. The AmnAB structure from the PDB (PDB ID: 8IHG) has an approximate radius of 4 nm, while the HO shell has a radius of about 20 nm. Under idealized conditions, where the enzymes are packed without any overlap and there was no solvent present, the shell could theoretically accommodate around 125 copies of AmnAB.† However, considering the macromolecular crowding constraint, where only about 20%⁶³ of the bacterial cytoplasmic volume is typically occupied by proteins, the realistic estimate would be approximately 25 copies of AmnAB within the shell interior. The reported k_cat for 2-AP conversion by AmnAB is 277 per s,¹⁹ suggesting that a fully loaded BMC could turn over slightly fewer than 7000 substrate molecules per second (approximately 1.6 × 10⁻²⁰ moles per s).

What conditions would support that kind of turnover? Clearly the inbound reactant flux and outbound product flux must match turnover, as otherwise the internal concentration would increase or decrease, changing the favorability for the reaction. Fick's law of diffusion relates molecular flux J to the concentration gradient between compartments (ΔC), dependent only on the surface area between compartments A and the permeability P.

J = PAΔC (4)

Eqn (4) can be rearranged to solve for a target concentration, . If we assume that the BMC shell is again approximately spherical with a internal radius of 20 nm, and that we take the minimum permeability we calculate from our simulations (0.1 cm s⁻¹), the concentration difference that would saturate encapsulated AmnAB is only 2.3 μM.‡

Is a concentration gradient in the μM range realistic? Compared to the 70 μM K_m for 2-AP to AmnAB,¹⁹ a concentration on the μM scale is small. Indeed, there are many metabolites where these concentrations are common, and for these metabolites our permeability values suggest that there is very little difference in the concentration inside and outside the shell. It could be that our assumptions are too generous for how many enzymes realistically fit inside a shell, but for HO shells like we are working with here, one really cannot be wrong by more than a factor of 10–20 without running into real physical limitations related to either creating a totally empty shell without any enzyme or an impossibly full shell.

The high predicted permeabilities are also consistent with the in vitro work in Fig. 7 and 8. While there are always questions for in vitro work if all the shells are always complete, the methods used to construct the shells in this study are quite robust.⁵⁴ Since the enzymes present in both systems are tethered within the shell, these enzymes are encapsulated and clearly can turn over while encapsulated. This demonstrates that the high permeabilities quantified through simulation are plausible.

With this thought experiment and our in vitro results in mind, our view is that BMC shells are not particularly selective in their permeability to the substrates in question. Instead of acting as a highly selective barrier, the BMC shell can improve catalytic efficiency by lengthening the path intermediates must travel before escaping, which increases their likelihood of being converted by encapsulated enzymes, as shown for CO₂ in carboxysomes.¹¹ This would reduce toxicity for intermediates present along the encapsulated reaction pathway by creating a kinetic barrier to transit.⁹ Table 1 and Fig. 5 show a narrow range of permeability coefficients for small molecules with varying charge and size, with less than a factor of 50 between the largest and smallest permeability.

That is not to say that varied permeability cannot be a contributor to BMC function. For instance, Park et al.¹⁵ computed a lower barrier to the passage of propanediol compared with the more reactive (and thus toxic) propionaldehyde, and estimated that the diol would permeate ten times faster than the aldehyde.¹⁵ However, our findings are much more muted, with a less than 2 fold selectivity for 2-AP over the reactive aldehyde, and the difference between the aldehyde and picolinate is even smaller (Table 1).

We suspect that the negative charges on most of our metabolites contribute to the high permeability. BMCs often have an asymmetric charge distributions, with one side often being more positive than the other. In CcmK1 and CcmK2 components of carboxysomes, the outer surface is mostly negative and the inner surface, including the pore, is positively charged.⁶⁴ This charge arrangement means that it is more favorable for the negatively charged photosynthetic metabolites to be inside rather than outside the shell, and likely helps attract and guide negatively charged metabolites through the pore.⁶⁴ The ARO system, which is involved in the native pathway to degrade 2-AP, can show oppositely patterned charge asymmetries, with either positive or negative charges predominating in the interior.¹⁹ The alternative charge arrangements in ARO are likely neutral to permeability in this case, as 2-AP itself is uncharged. The HO BMC used here as a model has a comparable charge imbalance to the carboxysome (Fig. S14), which may similarly facilitate transit for negatively charged substrates. Indeed, the positively charged ammonium cation is among the least permeable molecules in Table 1.

Conclusion

In this work we, have used atomistic molecular dynamics simulations to probe metabolite transports from two pathways catalyzed by distinct enzymes and encapsulated within a bacterial microcompartment (Fig. 1). The free energy, local diffusion and permeation coefficients were reported for each metabolite in Fig. 3 and 4. Our key finding is that the free energy barrier for most metabolites is comparable and small, suggesting no major restriction by the BMC while passing through the pores. Using the inhomogeneous solubility diffusion model, we find that the permeabilities for all metabolites are similar and rather high, comparable to a lipophilic molecule through a lipid bilayer. This indicates that BMC shell pores allow rapid metabolite passage with minimal selectivity. Compared to lipid membranes, diffusion is much faster, and similar to protein nanopores or viral capsid⁶⁵ channels. This efficient transport helps retain intermediates and supports catalytic efficiency, even without strict molecular selectivity. Thus, we anticipate no limitations to catalysis based on permeability, which is corroborated by fast catalysis for encapsulated enzymes in vitro. Instead, challenges to encapsulating enzymes within the shell itself may be the bigger challenge to developing designer microcompartments tailored to a specific reaction of interest.

Author contributions

NSY performed the computational modeling, molecular dynamics simulations, data analysis, and manuscript writing. YW and JFL carried out experiments to track activity for encapsulated AmnAB and NrfA, respectively. ELH, RPH, and JVV contributed to manuscript writing, review and editing, and funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

All input scripts to build and run molecular simulations are made publicly available on Zenodo (https://doi.org/10.5281/zenodo.16413720), along with analysis scripts. Full outputs are available upon request from the author, and are multiple tens of terabytes.

Supplementary information (SI): additional figures to support the manuscript, including extended discussion detailing how the small molecule parameters for the ligands were determined, free energy convergence estimation, and other ancilliary figures to the main manuscript. See DOI: https://doi.org/10.1039/d5nr03156g.

Acknowledgements

Research was supported as part of the Center for Catalysis in Biomimetic Confinement, an Energy Frontier Research Center funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), under Award DE-SC0023395. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, and additional computational support was provided by Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. This work was supported in part through computational resources and services provided by the Institute for Cyber-Enabled Research at Michigan State University.

References

1. M. Sutter, M. R. Melnicki, F. Schulz, T. Woyke and C. A. Kerfeld, Nat. Commun., 2021, 12, 3809.
2. A. M. Stewart, K. L. Stewart, T. O. Yeates and T. A. Bobik, Trends Biochem. Sci., 2021, 46, 406–416.
3. N. W. Kennedy, C. E. Mills, T. M. Nichols, C. H. Abrahamson and D. Tullman-Ercek, Curr. Opin. Microbiol., 2021, 63, 36–42.
4. S. M. Rose, A. Radhakrishnan and S. Sinha, J. Mater. Chem. B, 2023, 11, 4842–4854.
5. H. Kirst and C. A. Kerfeld, Biochem. Soc. Trans., 2021, 49, 1085–1098.
6. C. A. Kerfeld, C. Aussignargues, J. Zarzycki, F. Cai and M. Sutter, Nat. Rev. Microbiol., 2018, 16, 277–290.
7. E. J. Young, H. Kirst, M. E. Dwyer, J. V. Vermaas and C. A. Kerfeld, ACS Synth. Biol., 2025, 14, 1405–1413.
8. J. M. Ochoa and T. O. Yeates, Curr. Opin. Microbiol., 2021, 62, 51–60.
9. T. A. Bobik and A. M. Stewart, Curr. Opin. Microbiol., 2021, 62, 76–83.
10. B. J. Greber, M. Sutter and C. A. Kerfeld, Structure, 2019, 27, 749–763.
11. D. Sarkar, C. Maffeo, M. Sutter, A. Aksimentiev, C. A. Kerfeld and J. V. Vermaas, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2402277121.
12. M. Faulkner, I. Szabó, S. L. Weetman, F. Sicard, R. G. Huber, P. J. Bond, E. Rosta and L.-N. Liu, Sci. Rep., 2020, 10, 17501.
13. P. Mahinthichaichan, D. M. Morris, Y. Wang, G. J. Jensen and E. Tajkhorshid, J. Phys. Chem. B, 2018, 122, 9110–9118.
14. D. S. Trettel, C. Neale, M. Zhao, S. Gnanakaran and C. R. Gonzalez-Esquer, Sci. Rep., 2023, 13, 15738.
15. J. Park, S. Chun, T. A. Bobik, K. N. Houk and T. O. Yeates, J. Phys. Chem. B, 2017, 121, 8149–8154.
16. S. Raza, N. S. Yadav, A. Jussupow, C. A. Kerfeld, M. Feig and J. V. Vermaas, J. Phys. Chem. B, 2025, 129, 12811–12827.
17. J. S. Plegaria and C. A. Kerfeld, Curr. Opin. Biotechnol., 2018, 51, 1–7.
18. K. Hird, J. O. Campeciño, N. Lehnert and E. L. Hegg, J. Inorg. Biochem., 2024, 256, 112542.
19. L. Doron, M. Sutter and C. A. Kerfeld, mBio, 2023, e01216–e01223.
20. P. M. H. Kroneck, J. Biol. Inorg. Chem., 2022, 27, 1–21.
21. B. A. Donlon, E. Razo-Flores, G. Lettinga and J. A. Field, Biotechnol. Bioeng., 2000, 51, 439–449.
22. A. Turmo, C. R. Gonzalez-Esquer and C. A. Kerfeld, FEMS Microbiol. Lett., 2017, 364, fnx176.
23. C. S. Crowley, D. Cascio, M. R. Sawaya, J. S. Kopstein, T. A. Bobik and T. O. Yeates, J. Biol. Chem., 2010, 285, 37838–37846.
24. C. Chowdhury, S. Chun, A. Pang, M. R. Sawaya, S. Sinha, T. O. Yeates and T. A. Bobik, Proc. Natl. Acad. Sci. U. S. A., 2015, 112, 2990–2995.
25. H. Kirst and C. A. Kerfeld, BMC Biol., 2019, 17, 79.
26. A. Hagen, M. Sutter, N. Sloan and C. A. Kerfeld, Nat. Commun., 2018, 9, 2881.
27. R. M. Venable, A. Krämer and R. W. Pastor, Chem. Rev., 2019, 119, 5954–5997.
28. S. Raza, D. Sarkar, L. J. G. Chan, J. Mae, M. Sutter, C. J. Petzold, C. A. Kerfeld, C. Y. Ralston, S. Gupta and J. V. Vermaas, ACS Omega, 2024, 9, 35503–35514.
29. Y. Sugita, A. Kitao and Y. Okamoto, J. Chem. Phys., 2000, 113, 6042–6051.
30. C. T. Lee, J. Comer, C. Herndon, N. Leung, A. Pavlova, R. V. Swift, C. Tung, C. N. Rowley, R. E. Amaro, C. Chipot, Y. Wang and J. C. Gumbart, J. Chem. Inf. Model., 2016, 56, 721–733.
31. S. J. Marrink and H. J. C. Berendsen, J. Phys. Chem., 1996, 100, 16729–16738.
32. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
33. J. C. Phillips, D. J. Hardy, J. D. C. Maia, J. E. Stone, J. V. Ribeiro, R. C. Bernardi, R. Buch, G. Fiorin, J. Hénin, W. Jiang, R. McGreevy, M. C. R. Melo, B. K. Radak, R. D. Skeel, A. Singharoy, Y. Wang, B. Roux, A. Aksimentiev, Z. Luthey-Schulten, L. V. Kalé, K. Schulten, C. Chipot and E. Tajkhorshid, J. Chem. Phys., 2020, 153, 044130.
34. G. Fiorin, F. Marinelli and J. D. Faraldo-Gómez, J. Comput. Chem., 2020, 41, 449–459.
35. D. J. Sindhikara, D. J. Emerson and A. E. Roitberg, J. Chem. Theory Comput., 2010, 6, 2804–2808.
36. J. Huang, S. Rauscher, G. Nawrocki, T. Ran, M. Feig, B. L. De Groot, H. Grubmüller and A. D. MacKerell, Nat. Methods, 2017, 14, 71–73.
37. K. Vanommeslaeghe, E. Hatcher, C. Acharya, S. Kundu, S. Zhong, J. Shim, E. Darian, O. Guvench, P. Lopes, I. Vorobyov and A. D. Mackerell, J. Comput. Chem., 2010, 31, 671–690.
38. A. Kumar, O. Yoluk and A. D. MacKerell, J. Comput. Chem., 2020, 41, 958–970.
39. S. E. Feller, Y. Zhang, R. W. Pastor and B. R. Brooks, J. Chem. Phys., 1995, 103, 4613–4621.
40. S. Miyamoto and P. A. Kollman, J. Comput. Chem., 1992, 13, 952–962.
41. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
42. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
43. M. R. Shirts, D. L. Mobley, J. D. Chodera and V. S. Pande, J. Phys. Chem. B, 2007, 111, 13052–13063.
44. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
45. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, SciPy 1.0 Contributors, A. Vijaykumar, A. P. Bardelli, A. Rothberg, A. Hilboll, A. Kloeckner, A. Scopatz, A. Lee, A. Rokem, C. N. Woods, C. Fulton, C. Masson, C. Häggström, C. Fitzgerald, D. A. Nicholson, D. R. Hagen, D. V. Pasechnik, E. Olivetti, E. Martin, E. Wieser, F. Silva, F. Lenders, F. Wilhelm, G. Young, G. A. Price, G.-L. Ingold, G. E. Allen, G. R. Lee, H. Audren, I. Probst, J. P. Dietrich, J. Silterra, J. T. Webber, J. Slavič, J. Nothman, J. Buchner, J. Kulick, J. L. Schönberger, J. V. de Miranda Cardoso, J. Reimer, J. Harrington, J. L. C. Rodríguez, J. Nunez-Iglesias, J. Kuczynski, K. Tritz, M. Thoma, M. Newville, M. Kümmerer, M. Bolingbroke, M. Tartre, M. Pak, N. J. Smith, N. Nowaczyk, N. Shebanov, O. Pavlyk, P. A. Brodtkorb, P. Lee, R. T. McGibbon, R. Feldbauer, S. Lewis, S. Tygier, S. Sievert, S. Vigna, S. Peterson, S. More, T. Pudlik, T. Oshima, T. J. Pingel, T. P. Robitaille, T. Spura, T. R. Jones, T. Cera, T. Leslie, T. Zito, T. Krauss, U. Upadhyay, Y. O. Halchenko and Y. Vázquez-Baeza, Nat. Methods, 2020, 17, 261–272.
46. C. R. Harris, K. J. Millman, S. J. Van Der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. Van Kerkwijk, M. Brett, A. Haldane, J. F. Del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke and T. E. Oliphant, Nature, 2020, 585, 357–362.
47. A. L. Ferguson, J. Comput. Chem., 2017, 38, 1583–1605.
48. G. Hummer, New J. Phys., 2005, 7, 34.
49. J. V. Vermaas, G. J. Bentley, G. T. Beckham and M. F. Crowley, J. Phys. Chem. B, 2018, 122, 10349–10361.
50. J. V. Vermaas, G. T. Beckham and M. F. Crowley, J. Phys. Chem. B, 2017, 121, 11311–11324.
51. F. B. Sheinerman and C. L. Brooks, J. Mol. Biol., 1998, 278, 439–456.
52. J. V. Vermaas and E. Tajkhorshid, Biochemistry, 2017, 56, 281–293.
53. D. M. Boren, S. Kredi, E. Positselskaya, M. Giladi, Y. Haitin and J. V. Vermaas, Protein Sci., 2025, 34, e70167.
54. K. L. Range, T. K. Chiang, A. Pramanik, J. F. Landa, S. N. Snyder, X. Zuo, D. M. Tiede, L. M. Utschig, E. L. Hegg, M. Sutter, C. A. Kerfeld and C. Y. Ralston, ACS Nano, 2025, 19, 11913–11923.
55. J. Campeciño, S. Lagishetty, Z. Wawrzak, V. Sosa Alfaro, N. Lehnert, G. Reguera, J. Hu and E. L. Hegg, J. Biol. Chem., 2020, 295, 11455–11465.
56. M. Sutter, B. Greber, C. Aussignargues and C. A. Kerfeld, Science, 2017, 356, 1293–1297.
57. J. V. Vermaas, R. A. Dixon, F. Chen, S. D. Mansfield, W. Boerjan, J. Ralph, M. F. Crowley and G. T. Beckham, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 23117–23123.
58. S. Raza, M. Miller, B. Hamberger and J. V. Vermaas, J. Phys. Chem. B, 2023, 127, 1144–1157.
59. S. Raza, T. H. Sievertsen, S. Okumoto and J. V. Vermaas, Phytochemistry, 2024, 217, 113891.
60. S.-J. Marrink and H. J. C. Berendsen, J. Phys. Chem., 1994, 98, 4155–4168.
61. A. B. Mamonov, R. D. Coalson, M. L. Zeidel and J. C. Mathai, J. Gen. Physiol., 2007, 130, 111–116.
62. J. Tong, Z. Wu, M. M. Briggs, K. Schulten and T. J. McIntosh, Biophys. J., 2016, 111, 90–99.
63. M. Feig, R. Harada, T. Mori, I. Yu, K. Takahashi and Y. Sugita, J. Mol. Graphics Modell., 2015, 58, 1–9.
64. M. Sutter, S. McGuire, B. Ferlez and C. A. Kerfeld, ACS Synth. Biol., 2019, 8, 668–674.
65. C. Xu, D. K. Fischer, S. Rankovic, W. Li, R. A. Dick, B. Runge, R. Zadorozhnyi, J. Ahn, C. Aiken, T. Polenova, A. N. Engelman, Z. Ambrose, I. Rousso and J. R. Perilla, PLos Biol., 2020, 18, e3001015.

---

Footnotes

† .

‡ .

---