Advertisement

Molecular simulation of lignin-related aromatic compound permeation through gram-negative bacterial outer membranes

  • Josh V. Vermaas
    Correspondence
    For correspondence: Josh V. Vermaas; Michael F. Crowley; Gregg T. Beckham
    Affiliations
    Biosciences Center, National Renewable Energy Laboratory, Golden, Colorado, USA

    National Center for Computational Sciences, Oak Ridge National Laboratory, Oak Ridge, Tennessee, USA

    MSU-DOE Plant Research Laboratory, Michigan State University, East Lansing, Michigan, USA

    Department of Biochemistry and Molecular Biology, Michigan State University, East Lansing, Michigan, USA
    Search for articles by this author
  • Michael F. Crowley
    Correspondence
    For correspondence: Josh V. Vermaas; Michael F. Crowley; Gregg T. Beckham
    Affiliations
    Renewable Resources and Enabling Sciences Center, National Renewable Energy, Laboratory, Golden, Colorado, USA
    Search for articles by this author
  • Gregg T. Beckham
    Correspondence
    For correspondence: Josh V. Vermaas; Michael F. Crowley; Gregg T. Beckham
    Affiliations
    Renewable Resources and Enabling Sciences Center, National Renewable Energy, Laboratory, Golden, Colorado, USA
    Search for articles by this author
Open AccessPublished:October 20, 2022DOI:https://doi.org/10.1016/j.jbc.2022.102627
      Lignin, an abundant aromatic heteropolymer in secondary plant cell walls, is the single largest source of renewable aromatics in the biosphere. Leveraging this resource for renewable bioproducts through targeted microbial action depends on lignin fragment uptake by microbial hosts and subsequent enzymatic action to obtain the desired product. Recent computational work has emphasized that bacterial inner membranes are permeable to many aromatic compounds expected from lignin depolymerization processes. In this study, we expand on these findings through simulations for 42 lignin-related compounds across a gram-negative bacterial outer membrane model. Unbiased simulation trajectories indicate that spontaneous crossing for the full outer membrane is relatively rare at molecular simulation timescales, primarily due to preferential membrane partitioning and slow diffusion within the lipopolysaccharide layer within the outer membrane. Membrane partitioning and permeability coefficients were determined through replica exchange umbrella sampling simulations to overcome sampling limitations. We find that the glycosylated lipopolysaccharides found in the outer membrane increase the permeation barrier to many lignin-related compounds, particularly the most hydrophobic compounds. However, the effect is relatively modest; at industrially relevant concentrations, uncharged lignin-related compounds will readily diffuse across the outer membrane without the need for specific porins. Together, our results provide insight into the permeability of the bacterial outer membrane for assessing lignin fragment uptake and the future production of renewable bioproducts.

      Keywords

      Abbreviations:

      ABF (adaptive biasing force), IM (inner membrane), LPS (lipopolysaccharide), LRC (lignin related compound), OM (outer membrane), REUS (replica exchange umbrella sampling)
      Lignin, an aromatic heteropolymer found in terrestrial plant cell walls, is the largest single source of aromatics in the biosphere (
      • Boerjan W.
      • Ralph J.
      • Baucher M.
      Lignin biosynthesis.
      ). The substituted aromatics represent a high energy investment during biosynthesis, with substantial opportunities for valorization (
      • Ragauskas A.J.
      • Beckham G.T.
      • Biddy M.J.
      • Chandra R.
      • Chen F.
      • Davis M.F.
      • et al.
      Lignin valorization: improving lignin processing in the biorefinery.
      ). Utilizing these abundant lignin resources industrially to supplant fossil-derived products would be significant for sustainable economic development through a circular bioeconomy. In nature, fungal and bacterial enzymes breakdown the lignin polymer into smaller fragments (
      • Janusz G.
      • Pawlik A.
      • Sulej J.
      • Świderska-Burek U.
      • Jarosz-Wilkołazka A.
      • Paszczyński A.
      Lignin degradation: microorganisms, enzymes involved, genomes analysis and evolution.
      ,
      • Bugg T.D.H.
      • Ahmad M.
      • Hardiman E.M.
      • Rahmanpour R.
      Pathways for degradation of lignin in bacteria and fungi.
      ).
      Once fragmented, biological conversion for lignin and lignin related compounds (LRCs) through funneling into targeted end products has substantial promise to valorize lignin and other waste aromatics (
      • Borchert A.J.
      • Henson W.R.
      • Beckham G.T.
      Challenges and opportunities in biological funneling of heterogeneous and toxic substrates beyond lignin.
      ,
      • Linger J.G.
      • Vardon D.R.
      • Guarnieri M.T.
      • Karp E.M.
      • Hunsinger G.B.
      • Franden M.A.
      • et al.
      Lignin valorization through integrated biological funneling and chemical catalysis.
      ,
      • Becker J.
      • Wittmann C.
      From systems biology to metabolically engineered cells - an omics perspective on the development of industrial microbes.
      ,
      • Abdelaziz O.Y.
      • Brink D.P.
      • Prothmann J.
      • Ravi K.
      • Sun M.
      • Garćıa-Hidalgo J.
      • et al.
      Biological valorization of low molecular weight lignin.
      ,
      • Kamimura N.
      • Takahashi K.
      • Mori K.
      • Araki T.
      • Fujita M.
      • Higuchi Y.
      • et al.
      Bacterial catabolism of lignin-derived aromatics: new findings in a recent decade: update on bacterial lignin catabolism.
      ). Aromatic-catabolic bacteria, such as the soil bacterium Pseudomonas putida, have substantial metabolic flexibility to utilize an array of LRC molecules in their environment as their initial feedstock toward building designer products (
      • Becker J.
      • Wittmann C.
      From systems biology to metabolically engineered cells - an omics perspective on the development of industrial microbes.
      ,
      • Erickson E.
      • Bleem A.
      • Kuatsjah E.
      • Werner A.Z.
      • DuBois J.L.
      • McGeehan J.E.
      • et al.
      Critical enzyme reactions in aromatic catabolism for microbial lignin conversion.
      ,
      • Beckham G.T.
      • Johnson C.W.
      • Karp E.M.
      • Salvachúa D.
      • Vardon D.R.
      Opportunities and challenges in biological lignin valorization.
      ,
      • Johnson C.W.
      • Salvachúa D.
      • Rorrer N.A.
      • Black B.A.
      • Vardon D.R.
      • St. John P.C.
      • et al.
      Innovative chemicals and materials from bacterial aromatic catabolic pathways.
      ). This flexibility could be enhanced further by integrating fungal pathways for lignin catabolism into microbial systems (
      • del Cerro C.
      • Erickson E.
      • Dong T.
      • Wong A.R.
      • Eder E.K.
      • Purvine S.O.
      • et al.
      Intracellular pathways for lignin catabolism in white-rot fungi.
      ).
      One underexplored aspect to support this work is how LRCs permeate into bacterial cells. P. putida and other gram-negative bacteria have two membranes, an inner membrane (IM) composed of phospholipids and an outer membrane (OM) that features a lipopolysaccharide (LPS) membrane leaflet (
      • Silhavy T.J.
      • Kahne D.
      • Walker S.
      The bacterial cell envelope.
      ). For these bacteria to catabolize LRCs, the compounds must traverse both of these membranes, either passively or via transport protein, and also diffuse through the peptidoglycan cell wall.
      Past computational work has demonstrated that uncharged LRCs passively permeate through a bacterial IM at rates commensurate with cellular metabolism (
      • Vermaas J.V.
      • Dixon R.A.
      • Chen F.
      • Mansfield S.D.
      • Boerjan W.
      • Ralph J.
      • et al.
      Passive membrane transport of lignin-related compounds.
      ). Additionally, vesicles that bud from the OM have been observed to catalyze lignin catabolism, suggesting that LRCs permeate the OM through an unresolved mechanism (
      • Salvachúa D.
      • Werner A.Z.
      • Pardo I.
      • Michalska M.
      • Black B.A.
      • Donohoe B.S.
      • et al.
      Outer membrane vesicles catabolize lignin-derived aromatic compounds in Pseudomonas putida KT2440.
      ). The OM LPS glycosylations within the outer leaflet form a hydrophilic region that could act as a permeation barrier to largely hydrophobic LRCs. If the permeation barrier is significant, transport proteins would be required to facilitate LRC flux across the OM. Alternatively, LRCs may passively permeate across the OM, just as was predicted for the IM (
      • Vermaas J.V.
      • Dixon R.A.
      • Chen F.
      • Mansfield S.D.
      • Boerjan W.
      • Ralph J.
      • et al.
      Passive membrane transport of lignin-related compounds.
      ) and recently validated experimentally for plant-like membranes (
      • Perkins M.L.
      • Schuetz M.
      • Unda F.
      • Chen K.T.
      • Bally M.B.
      • Kulkarni J.A.
      • et al.
      Monolignol export by diffusion down a polymerization-induced concentration gradient.
      ).
      Through molecular simulation tools, we test the passive permeation hypothesis across the OM by quantifying the permeability for the OM of P. putida. Molecular dynamics (MD) is a useful tool for studying small molecule interactions with lipid bilayers (
      • Martinotti C.
      • Ruiz-Perez L.
      • Deplazes E.
      • Mancera R.L.
      Molecular dynamics simulation of small molecules interacting with biological membranes.
      ) and comes with well-developed theory for determining permeability in silico (
      • Venable R.M.
      • Krämer A.
      • Pastor R.W.
      Molecular dynamics simulations of membrane permeability.
      ). While software for building LPS models has only recently become available (
      • Jo S.
      • Wu E.L.
      • Stuhlsatz D.
      • Klauda J.B.
      • MacKerell A.D.
      • Widmalm G.
      • et al.
      Lipopolysaccharide membrane building and simulation.
      ), molecular simulation has quickly proven to be an effective tool for studying LPS-bearing membranes, providing a detailed view into the molecular interactions within LPS and its environment (
      • Khalid S.
      • Piggot T.J.
      • Samsudin F.
      Atomistic and coarse grain simulations of the cell envelope of gram-negative bacteria: what have we learned?.
      ,
      • Ĺopez C.A.
      • Zgurskaya H.
      • Gnanakaran S.
      Molecular characterization of the outer membrane of Pseudomonas aeruginosa.
      ).
      Through equilibrium and nonequilibrium MD simulations for selected LRC compounds (Fig. 1), we determine that the OM can be a bigger barrier to permeation than the IM, particularly for the most lipophilic LRCs tested. However, the barrier is not sufficiently large to require protein-specific transport processes. The OM slows permeation by a few orders of magnitude compared to IM permeation, with the resulting permeability coefficients remaining high enough to support robust flux across the OM. The LPS region of the OM was found to hinder permeation the most, as the slow LPS dynamics reduce diffusion in that region significantly. Taken together, the results suggest that specific transport proteins are not required for most LRCs to permeate the OM within a biorefinery context, as potential passive fluxes in these scenarios exceed cellular LRC catabolic capacity.
      Figure thumbnail gr1
      Figure 1Lignin-related compounds used in this study. The aromatic molecules tested represent monomeric lignin compounds with diverse lignin functionalities. Each colored box groups together molecule functionality at position 1 of the aromatic ring. The box color is consistent with the color assigned to these functional groups in other Figures. Within the colored boxes, the lignin compounds vary based on hydroxylation and methoxylation patterns for H-, G-, S-, and C-type lignin molecules.

      Results

      Equilibrium analysis

      The unbiased equilibrium trajectories can provide insight into the typical behavior for small molecules near our OM model that is simultaneously glycosylated and lipidated. While individual small molecules can have radically different traces based on the stochastic processes underlying membrane insertion and translocation (Figs. S1–S41), across all 42 compounds under study, we observe only a single instance where a small molecule passively permeates across the OM model generated here within our 400 ns simulation duration. For the syringol system, a single syringol molecule (purple line, Fig. 2) is observed to enter into the glycosylated region marking the OM boundary. Approximately 150 ns into the simulation, this single syringol molecule enters more deeply into the glycosylated region, eventually reaching the lipid core. The syringol molecule remains in the outer leaflet for approximately 50 ns prior to a rapid exchange to the inner leaflet, where it remains for approximately 100 ns before briefly entering into solution and completing a full membrane transit event (Supplementary Animation 1).
      Figure thumbnail gr2
      Figure 2Traces for the 10 syringol molecules placed into a membrane context. Each uniquely colored traces mark the pathway for a single syringol molecule, drawn to correct for trajectories that cross the periodic boundary. To assist the reader in tracking the motion across the different compartments, the hydrophobic lipid core region of the membrane has a gray background, while the glycosylated region of the membrane has a yellow background color. Regions of aqueous solution have no background coloration. The purple trace is observed to permeate from the top of the glycolipid leaflet and cross the hydrophobic membrane thickness, briefly even entering solution on the far side.
      The rarity for full transits in our dataset supports the hypothesis that OMs pose a larger barrier to permeation than IMs for some compounds. In prior IM studies, phenol, guaiacol, and syringol readily translocated between leaflets in 500 ns equilibrium simulations (
      • Vermaas J.V.
      • Dixon R.A.
      • Chen F.
      • Mansfield S.D.
      • Boerjan W.
      • Ralph J.
      • et al.
      Passive membrane transport of lignin-related compounds.
      ), with multiple translocation events observed. Leaflet exchanges are also observed in the current membrane models, with the same three molecules readily crossing from the inner, unglycosylated leaflet into the opposing leaflet and back again (Figs. 2, S1, and S2). The missing component required for permeation are crossing events across the glycosylated region within the LPS-filled outer leaflet. While our dataset only contains a single permeation event, its existence within 100 microseconds of aggregated equilibrium sampling across all compounds suggests that the permeation rate across the OM may be fast on a biological timescale.
      Mechanistically, slow permeation across the LPS leaflet may be related to the slow diffusion processes observed in that leaflet during the equilibrium trajectories. Animation S1 already hints that the dynamics within LPS leaflet are slower than in the phospholipid leaflet. Water diffusion parallel and perpendicular to the membrane normal axis is slower within the LPS leaflet than within the phospholipid leaflet (Fig. 3). Slower water diffusion within the LPS leaflet is not limited to the LPS glycosylations but also extends to membrane depths where LPS acyl tails predominate. Mechanistically, the dense glycosylations may order and trap water or potentially other similar small molecules. In a different measure of lipid dynamics, the LPS leaflet in our OM model diffuses nearly 100 times more slowly than the phospholipid leaflet (Table 1), in line with prior computational estimates for the independent OM leaflets (
      • Piggot T.J.
      • Holdbrook D.A.
      • Khalid S.
      Electroporation of the E. coli and S. Aureus membranes: molecular dynamics simulations of complex bacterial membranes.
      ,
      • Soares T.
      • Straatsma T.
      Assessment of the convergence of molecular dynamics simulations of lipopolysaccharide membranes.
      ,
      • Balusek C.
      • Gumbart J.C.
      Role of the native outer-membrane environment on the transporter BtuB.
      ). The slow lateral lipid diffusion generally is caused by the increased mass for LPS molecules and the abundant polar interactions within the LPS glycosylations. The slow diffusion processes caused by abundant interactions up and down the LPS chain slow all dynamics in the LPS leaflet, retarding membrane crossing.
      Figure thumbnail gr3
      Figure 3Position-dependent water diffusion coefficients both parallel to (Waterk) and perpendicular to (Water) the membrane normal axis. For context, the black curves indicating the diffusion coefficient are overlaid on a molecular graphics representation for the glycosylated membrane following the color scheme from .
      Table 1Lateral diffusion coefficients for LPS and phospholipid leaflets in our OM model
      LeafletDlateral
      LPS5.6 ± 0.3 × 10−10 cm2s−1
      Phospholipid4.9 ± 0.2 10−8cm2s−1
      Notably, the slow dynamics in the LPS leaflet does not imply that the LPS leaflet is more ordered. Comparing between LPS and phospholipid tails within our model membrane indicates that the LPS acyl tails have a lower average order parameter at the equivalent carbon positions than their phospholipid counterparts (Fig. 4), indicating greater conformational variability in the LPS acyl tail than in the phospholipid tails. While prior LPS simulations did not explicitly compare acyl tail ordering with phospholipid acyl tails (
      • Luna E.
      • Kim S.
      • Gao Y.
      • Widmalm G.
      • Im W.
      Influences of Vibrio cholerae lipid A types on LPS bilayer properties.
      ,
      • Wu E.L.
      • Cheng X.
      • Jo S.
      • Rui H.
      • Song K.C.
      • Dávila-Contreras E.M.
      • et al.
      CHARMM-GUI membrane builder toward realistic biological membrane simulations.
      ), the magnitude and trends for -SCD are consistent with prior results both for LPS (
      • Luna E.
      • Kim S.
      • Gao Y.
      • Widmalm G.
      • Im W.
      Influences of Vibrio cholerae lipid A types on LPS bilayer properties.
      ,
      • Wu E.L.
      • Cheng X.
      • Jo S.
      • Rui H.
      • Song K.C.
      • Dávila-Contreras E.M.
      • et al.
      CHARMM-GUI membrane builder toward realistic biological membrane simulations.
      ) and phospholipids (
      • Piggot T.J.
      • Allison J.R.
      • Sessions R.B.
      • Essex J.W.
      On the calculation of acyl chain order parameters from lipid simulations.
      ). As we are averaging over multiple different acyl chains with a mixed lipid composition, the order parameter reduction near common unsaturation sites is less evident in Figure 4 than is typical from homogeneous lipid compositions that study each acyl tail independently.
      Figure thumbnail gr4
      Figure 4Mean acyl tail order parameter (-SCD) computed for all acyl tails in the LPS (red) and phospholipid (black) leaflet. SCD=<3cos2(θCH)1>/2 where θCH is the angle between a C-H bond vector from the trajectory and the membrane-normal axis (
      • Klauda J.B.
      • Venable R.M.
      • Freites J.A.
      • O’Connor J.W.
      • Tobias D.J.
      • Mondragon-Ramirez C.
      • et al.
      Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types.
      ). As the number of samples behind each datapoint is so high across all 16.8μs of equilibrium simulation, the SEM for each acyl tail order parameter does not rise outside the drawn point. LPS, lipopolysaccharide.

      Predicting permeation coefficients

      While equilibrium trajectory analysis serves to provide mechanistic insight into LRC translocation through bacterial OMs, the equilibrium trajectories alone are insufficient to compute accurate free energy profiles or kinetics for OM crossing. Without an applied biasing potential, the LRCs extensively and unproductively sample the reaction coordinate near membrane lipid tails where their interactions are most favorable. Through biased simulation, we are able to sample evenly across the bilayer translocation reaction coordinate, collecting statistics on rarely visited transition regions critical to determining the kinetics for OM permeation events.
      In the inhomogeneous solubility diffusion model (Equation 3) (
      • Marrink S.-J.
      • Berendsen H.J.C.
      Simulation of water transport through a lipid membrane.
      ,
      • Lee C.T.
      • Comer J.
      • Herndon C.
      • Leung N.
      • Pavlova A.
      • Swift R.V.
      • et al.
      Simulation-based approaches for determining membrane permeability of small compounds.
      ), the permeability coefficient depends only on the free energy and diffusivity profiles across a region of interest. The two essential profiles for selected G-type lignin monomers from our larger LRC set are provided in Figure 5 and indicate substantial changes in permeability depending on the LRC. The highest free energy barrier overall is observed for vanillate near the membrane midplane, reflecting the energetic cost for dehydrating the vanillate anion. If the carboxylate anion were to protonate, as in vanillic acid, the energetic penalty for crossing the membrane is substantially reduced. For many of the other LRCs tested, the largest barrier to permeation was observed at the interface between the lipid and polysaccharide portions for the LPS leaflet. At this interface, typical free energies for uncharged species were between 2-6 kcal mol−1 relative to solution (Figs. 5 and S42–S51), representing a significant but not insurmountable barrier to permeation.
      Figure thumbnail gr5
      Figure 5Free energy (top) and diffusivity (bottom) profiles for selected G-type lignin monomers included in our test set (). Each compound has two lines associated with the compound, a solid line reflecting the instantaneous best estimate for the quantity of interest and a dashed line indicating the spline fit used to numerically integrate Equation  to tabulate permeability coefficients in . For context, the plots are underlaid with a molecular representation for the glycosylated membrane, following the color scheme from . Free energy and diffusivity profiles for other LRCs are provided in .
      The free energy for uncharged LRCs was typically negative within the membrane core relative to solution, suggesting favorable partitioning for these small molecules into the membrane and rapid, low-barrier exchanges between leaflets. The free energy profiles also suggest that the preferred insertion depth for an LRC was leaflet dependent, with the LPS leaflet demonstrating a preferred insertion depth nearer to the membrane midplane than the phospholipid leaflet. The shorter acyl tail lengths and lower lipid ordering for LPS membrane components means that carbonyl-adjacent regions for the membrane are closer to the midplane for an LPS-containing leaflet than for our phospholipid model. Likewise, other equivalent condition regions within the membrane core exist and are shifted with respect to the membrane midplane. Another common feature for the LRC free energy profiles are barriers near the transition between the LPS and phospholipid leaflets. We ascribe the offset between the membrane center and the free energy maxima to the dynamical differences between the LPS and phospholipid leaflets, as the membrane center was defined based on the center of mass for the acyl tail termini.
      The diffusion coefficient profiles are similarly asymmetric in a leaflet-dependent manner. The diffusivity along the membrane normal within the phospholipid leaflet compares favorably with prior measures for diffusion coefficients along membrane transfer axes (
      • Vermaas J.V.
      • Dixon R.A.
      • Chen F.
      • Mansfield S.D.
      • Boerjan W.
      • Ralph J.
      • et al.
      Passive membrane transport of lignin-related compounds.
      ). For phospholipid membranes, the aqueous diffusion coefficients for LRCs are near 80Å2/ns, exhibiting a minima near the membrane/water interface before increasing LRC diffusion near the membrane center. The diffusive behavior within phospholipid leaflets, including diffusion minima near interfaces, have been noted previously for gas diffusion within membranes (
      • Ghysels A.
      • Venable R.M.
      • Pastor R.W.
      • Hummer G.
      Position-dependent diffusion tensors in anisotropic media from simulation: oxygen transport in and through membranes.
      ,
      • De Vos O.
      • Venable R.M.
      • Van Hecke T.
      • Hummer G.
      • Pastor R.W.
      • Ghysels A.
      Membrane permeability: characteristic times and lengths for oxygen and a simulation-based test of the inhomogeneous solubility-diffusion model.
      ), although alternative analysis methods yield alternative diffusivity profiles (
      • Lee C.T.
      • Comer J.
      • Herndon C.
      • Leung N.
      • Pavlova A.
      • Swift R.V.
      • et al.
      Simulation-based approaches for determining membrane permeability of small compounds.
      ,
      • Sicard F.
      • Koskin V.
      • Annibale A.
      • Rosta E.
      Position-dependent diffusion from biased simulations and markov state model analysis.
      ). Qualitatively, the trends are similar for the LPS leaflet. In the LPS leaflet, the diffusion coefficient is high in aqueous solution and somewhat lower at the membrane center. The slowest diffusion occurs near the membrane-glycosylation interface as the LRCs transition between lipophilic and hydrophobic regions. The diffusion at the membrane-glycosylation interface is approximately half the rate as what was observed in the membrane-water interface for the phospholipid leaflet (Fig. 5). The diffusive reduction in LPS is comparable in magnitude to the water diffusion parallel to the membrane normal axis presented in Figure 3.
      With the free energy and diffusivity profiles in hand, we are able to compute partitioning and permeability coefficients for the diverse LRCs in our simulation set (Tables 2 and S1) based on Equations 3 and 4. The partition coefficients are nearly uniformly positive, indicating that the LRCs partition into the membrane, with only a few charged molecules preferring solution (Table S1). The permeability coefficients vary by up to 13 orders of magnitude, with charged compounds such as vanillate demonstrating the slowest permeation. For uncharged species, particularly for common compounds enumerated in Table 2, permeability was often limited by permeation through the glycosylation attached to the LPS leaflet (Pmg), consistent with the identified regions of high free energy and low diffusion (Fig. 5).
      Table 2Partition (P) and Permeability (Pm) coefficients for G-type LRCs shown in Figure 5 in the OM mimic
      Compound namelogPlogPm(cms−1)logPmu (cms−1)logPmc (cms−1)logPmg (cms−1)
      Guaiacol4.0-0.8-1.9-0.3-4.8
      Vanillyl alcohol1.0-0.40.8-0.9-1.3
      Vanillin2.5-1.6-0.5-1.4-4.1
      Coniferyl alcohol1.9-1.0-0.0-0.4-2.9
      Vanillate1.3-10.7-0.9-12.01.0
      The permeability coefficient of crossing the entire membrane, going from aqueous solution to aqueous solution, is decomposed into a crossing permeability (Pmc) and extraction permeabilities into solution through the glycosylated (Pmg) and unglycosylated (Pmu) sides. The decomposition is done by adjusting the integrated bounds in Equation 3 to integrate between free energy minima found near the lipid-water and lipid-glycosylation interfaces. This is described by Equation 5 in the methods. The same information for all tested compounds is available in Table S1. Typical uncertainties for the tabulated quantities are 0.2 log units.

      Discussion

      Based exclusively on the results presented, the initial hypothesis that the LPS layer in a P. putida OM reduces LRC permeability appears plausible. Spontaneous permeation is rare within our simulation set (Fig. 2), and the glycosylations frequently present the largest barrier to permeation (Fig. 5 and Table 2). To demonstrate what practical impact the OM might have on LRC uptake on P. putida, we take the opportunity here to place the OM permeation coefficients into their biological context.

      Double membrane permeability and LRC chemistry

      Our model for LRC uptake in P. putida is quite simple, consisting of a LRC source external to the cell derived from biomass breakdown, an LRC sink representing internal metabolism, and the P. putida OM and IM (Fig. 6). The modeled passive transport system reaches a steady state where the molecular flux into the cell matches consumption by cellular metabolism. At steady state, Jo→m = Jm→i = Jo→i, where the effective permeability coefficient for the total process (Pmt) can be calculated from the permeability across the OM (Pmo) and IM (Pmi).
      Pmt=PmoPmiPmo+Pmi
      (1)


      Figure thumbnail gr6
      Figure 6Schematic for the simplified membrane transport processes through both the OM and IM in P. putida, highlighting the three primary aqueous compartments featuring LRC concentrations. Biomass breakdown increases LRC concentration outside the cell ([LRCo]), while aromatic catabolic processes inside the cell reduce the LRC concentration ([LRCi]). This establishes a concentration gradient that approaches steady state when the fluxes across the individual membranes (Jo→m and Jm→i). In this schematic, the hydrophobic membrane regions are represented by pink rectangles, while the OM LPS glycosylations are represented by yellow hexagons. IM, inner membrane; OM, outer membrane; LRC, lignin related compound; LPS, lipopolysaccharide.
      Based on the structure of Equation 1, Pmt is approximately half the permeation coefficient for either membrane when Pmo and Pmi are similar in magnitude. Otherwise, Pmt will be limited by the smaller of Pmo and Pmi.
      Thus, when looking at the larger passive transport process, the results from Tables 2 and S1 need to be compared with prior results for the IM (
      • Vermaas J.V.
      • Dixon R.A.
      • Chen F.
      • Mansfield S.D.
      • Boerjan W.
      • Ralph J.
      • et al.
      Passive membrane transport of lignin-related compounds.
      ). From Figure 7, we observe variation in relative permeation coefficients dependent on LRC chemistries. LRCs with additional oxygenation sites, such as cinnamic acids, are limited by the IM permeability rather than the OM permeability, which can be up to 1000 times larger. More lipophilic LRCs may have IM permeability coefficients 100 times larger than the OM permeability coefficient. These lipophilic LRCs include phenols found in aqueous waste streams from catalytic fast pyrolysis (
      • Wilson A.N.
      • Dutta A.
      • Black B.A.
      • Mukarakate C.
      • Magrini K.
      • Schaidle J.A.
      • et al.
      Valorization of aqueous waste streams from thermochemical biorefineries.
      ). Other compounds have computed permeability coefficients that are comparable in magnitude. Thus, even though the largest barrier to crossing the OM is typically within the LPS glycosylation region, the overall impact on passive permeation coefficient across both bilayers is relatively low. The permeation barrier provided by the LPS glycosylation may reduce Pmt by 3 orders of magnitude. By comparison, changing LRC chemistry to include a negative charge reduces Pmt by approximately 10 orders of magnitude, a far more significant change with more far-reaching biological consequences.
      Figure thumbnail gr7
      Figure 7Outer-versus inner-membrane permeability for compounds within our test set () where IM permeability had been previously calculated (
      • Vermaas J.V.
      • Dixon R.A.
      • Chen F.
      • Mansfield S.D.
      • Boerjan W.
      • Ralph J.
      • et al.
      Passive membrane transport of lignin-related compounds.
      ). To guide the eye, we have added a dashed line to indicate where OM Pm = IM Pm. Below this line, Pmo < Pmi, indicating that the OM is rate limiting. The converse is true above the line. IM, inner membrane; OM, outer membrane.
      By regrouping the information from Figure 7 by the LRC type from Figure 1, we can evaluate how hydroxylation and methoxylation for LRC compounds change permeability (Fig. 8). From Figure 8, we find that methoxylated S- and G-type LRCs have distributions shifted to the left in Figure 8 compared with H- and C-type LRCs that are only hydroxylated. The shifted distribution indicates that the methoxy-substituted LRCs are less permeable across the OM than their hydroxylated counterparts. Hydroxyl groups can act as hydrogen bonding donors and acceptors, unlike methoxy groups, which can only act as hydrogen bonding acceptors. Within the LPS region, we find that hydroxyl groups form significantly more hydrogen bonds per functional group (Fig. S52). The favorable interactions with the LPS region lower the permeation barrier for hydroxylated LRCs. Thus, one could hypothesize that hydroxylated LRCs may be catabolized more quickly by P. putida.
      Figure thumbnail gr8
      Figure 8Replotting the data in by the type of lignin monomer as a cumulative distribution function, determining whether the small molecules are more likely to permeate through the OM or the IM. IM, inner membrane; OM, outer membrane; LRC, lignin related compound; LPS, lipopolysaccharide.

      Practical flux calculations

      With a permeability coefficient in hand for the combined OM and IM system based on Equation 1, we can calculate potential fluxes based on the permeability coefficient, the area for the permeating surface A, and the concentration gradient across the dividing surface.
      Joi=PmtA([LRC]o[LRC]i)
      (2)


      The key surface area and concentration gradient parameters within Equation 2 can be estimated from experimental sources. High mM titers for phenolic LRCs have been observed in prior microbial engineering efforts to produce phenols (
      • Wynands B.
      • Lenzen C.
      • Otto M.
      • Koch F.
      • Blank L.M.
      • Wierckx N.
      Metabolic engineering of Pseudomonas taiwanensis VLB120 with minimal genomic modifications for high-yield phenol production.
      ,
      • Wierckx N.J.P.
      • Ballerstedt H.
      • de Bont J.A.M.
      • Wery J.
      Engineering of solvent-tolerant Pseudomonas putida S12 for bioproduction of phenol from glucose.
      ) and would represent an upper bound for the concentration difference within Equation 2 P. putida could tolerate. The membrane surface area for P. putida can be estimated from micrographs showing a rod-like shape with long axis of approximately 1 μm and a 0.4 μm diameter (
      • Salvachúa D.
      • Werner A.Z.
      • Pardo I.
      • Michalska M.
      • Black B.A.
      • Donohoe B.S.
      • et al.
      Outer membrane vesicles catabolize lignin-derived aromatic compounds in Pseudomonas putida KT2440.
      ), yielding 1.25 μm2. Assuming a 1 mM concentration difference, the net inward flux for phenol predicted by Equation 2 is 25 fmol/s (1.5 × 1010 molecules/s), limited by Pmo=2 cm/s. While the OM limits flux and reduces permeation 10-fold compared with prior IM simulations (
      • Vermaas J.V.
      • Dixon R.A.
      • Chen F.
      • Mansfield S.D.
      • Boerjan W.
      • Ralph J.
      • et al.
      Passive membrane transport of lignin-related compounds.
      ), the potential permeation rate remains very high in this example.
      Many LRCs tested will exhibit similar behavior, as their permeabilities are broadly similar to that of phenol. Among the uncharged LRCs, sinapyl alcohol has the smallest permeability coefficient, with Pm= 10−2.9 cm/s, around 1600 times less permeable than the phenol example above, yielding permeation rates on the order of 107 molecules/s. With this high potential flux, the actual flux is likely controlled by the steady state rate of LRC synthesis or catabolism. Enzymes involved in aromatic demethylation turn over few times a second (
      • Mallinson S.J.B.
      • Machovina M.M.
      • Silveira R.L.
      • Garcia-Borra`s M.
      • Gallup N.
      • Johnson C.W.
      • et al.
      A promiscuous cytochrome P450 aromatic O-demethylase for lignin bioconversion.
      ,
      • Notonier S.
      • Werner A.Z.
      • Kuatsjah E.
      • Dumalo L.
      • Abraham P.E.
      • Hatmaker E.A.
      • et al.
      Metabolism of syringyl lignin-derived compounds in Pseudomonas putida enables convergent production of 2-pyrone-4,6-dicarboxylic acid.
      ), and lignin acetylation enzymes turn over in tens of seconds (
      • de Vries L.
      • MacKay H.A.
      • Smith R.A.
      • Mottiar Y.
      • Karlen S.D.
      • Unda F.
      • et al.
      HBMT1, a BAHD-family monolignol acyltransferase, mediates lignin acylation in poplar.
      ). Slow lignin metabolism implies that passive flux across the OM is not rate limiting, even in the absence of specific transport proteins. Instead, the concentration gradient implied by Equation 2 for uncharged LRCs would shrink such that the flux matches metabolic processes.
      For charged LRCs, permeability is substantially slower. If we replace phenol from our example above with p-hydroxy-benzoate, the permeation rate would slow to a few molecules per hour. For such a small flux, ingress and egress control by transport proteins is highly likely, tempered by the possibility that many permeation events for carboxylates may occur when the carboxylate protonates to form the neutral acid. The acids demonstrate permeability coefficients that are significantly faster and thus may represent the major species that passively permeates, particularly in acidic conditions.

      Model limitations

      While the main message that OM permeability is actually quite fast for LRCs is clear, there are potential areas of concern for our models. The first is how well the inhomogeneous solubility diffusion model and the underlying simulation methods and force fields recapitulate reality. Indeed, systematic surveys suggest that molecular simulation can overestimate permeability by an order of magnitude (
      • Lee C.T.
      • Comer J.
      • Herndon C.
      • Leung N.
      • Pavlova A.
      • Swift R.V.
      • et al.
      Simulation-based approaches for determining membrane permeability of small compounds.
      ). In our view, a modest systematic bias towards accelerated permeability is acceptable, as the potential passive fluxes are sufficiently large that overestimating these by a factor of 10 does not materially change whether permeation or metabolism limits the steady state flux through the system.
      However, a larger area of uncertainty arises from our choice for a relatively thin LPS layer. Removing the O-antigen entirely was necessary from a computational perspective to keep the simulation sizes tractable for a large set of LRCs. Consequently, we may overestimate permeability through the LPS layer, as additional glycosylation may offer greater resistance to permeation. If we make the assumption that the free energy and diffusivity across the additional O-antigen thickness can be taken from point where the free energy reaches a maximum within the LPS (Equation 6), we find that the O-antigen thickness required to reduce permeability tenfold is between 2.8 and 14.2 nm for uncharged LRCs. Since additional thickness would only linearly add resistance to permeation (
      • Marrink S.-J.
      • Berendsen H.J.C.
      Simulation of water transport through a lipid membrane.
      ,
      • Lee C.T.
      • Comer J.
      • Herndon C.
      • Leung N.
      • Pavlova A.
      • Swift R.V.
      • et al.
      Simulation-based approaches for determining membrane permeability of small compounds.
      ), it is likely infeasible for realistic O-antigen lengths to reduce OM permeability by more than 10 to 100 times unless we underestimate the free energy barrier in our current model construction.
      To test this hypothesis, we directly compare free energy and diffusivity profiles determined through adaptive biasing force (ABF) calculations for guaiacol in a membrane where the O-antigen is present to the original replica exchange umbrella sampling (REUS) simulations that lacked the O-antigen (Fig. 9). Overall, the free energy profile has a similar structure, with minima near the acyl tail carbonyls and a maxima where the LPS leaflet transitions from a lipophilic to a hydrophilic environment. Crucially, the barrier within the free energy profile has a similar maximum and returns to zero in the glycosylation region both with and without O-antigen present. Prior studies that calculate free energy barriers across LPS membranes arrive at similar small molecule free energy profiles (
      • Carpenter T.S.
      • Parkin J.
      • Khalid S.
      The free energy of small solute permeation through the Escherichia coli outer membrane has a distinctly asymmetric profile.
      ), suggesting that indeed the interface between the hydrophobic and hydrophilic regions within LPS are the primary barrier to permeation for nonpolar compounds. The extended sampling within the ABF calculation appears to have broadly improved the symmetry for the profile within the membrane.
      Figure thumbnail gr9
      Figure 9Free energy and diffusivity profile comparison for guaiacol. Two profiles are shown here, one where the membrane lacked an O-antigen and the profiles were computed via REUS calculations (black line, used throughout the text) and another profile set computed with a membrane containing an O-antigen, where the profile was computed via ABF (red line). Dashed lines represent the spline fits used to integrate Equation . ABF, adaptive biasing force; REUS, replica exchange umbrella sampling.
      The diffusivity profiles are substantially different between the two methodologies, which have been noted previously when comparing REUS and ABF results (
      • Lee C.T.
      • Comer J.
      • Herndon C.
      • Leung N.
      • Pavlova A.
      • Swift R.V.
      • et al.
      Simulation-based approaches for determining membrane permeability of small compounds.
      ). We also note that our ABF sampling strategy created discontinuities within the diffusivity profile for the O-antigen membrane, which we attempt to correct via a spline fit (Fig. 9). In truth, the diffusivity matters very little in the overall permeability picture. Remembering that the computed guaiacol permeability coefficient is 10−0.76 cm/s when the O-antigen is absent (rounds to −0.8 in Table 2), we integrate Equation 3 to estimate permeability when the O-antigen is present. The equivalent coefficients computed from the ABF results are either 10−0.81 cm/s (using REUS diffusivity and ABF free energy) or 10−0.40 cm/s (using ABF for both diffusivity and free energy). Based on this limited evidence, adding or removing the O-antigen on LPS does not materially change the permeability for the OM as a whole for LRCs.
      While the permeability through the OM and IM for uncharged LRCs is significant, the inherent assumption behind Equation 2 is that the peptidoglycan layer sandwiched between the OM and IM contributes little resistance to LRC permeability. While the peptidoglycan layer in P. putida is in an aqueous phase, interactions between the dense peptidoglycan mesh (
      • Gumbart J.C.
      • Beeby M.
      • Jensen G.J.
      • Roux B.
      Escherichia coli peptidoglycan structure and mechanics as predicted by atomic-scale simulations.
      ) and the LRCs may impede permeation. Strong glycan–lignin interactions have been observed in prior simulations (
      • Vermaas J.V.
      • Petridis L.
      • Qi X.
      • Schulz R.
      • Lindner B.
      • Smith J.C.
      Mechanism of lignin inhibition of enzymatic biomass deconstruction.
      ,
      • Vermaas J.V.
      • Crowley M.F.
      • Beckham G.T.
      A quantitative molecular atlas for interactions between lignin and cellulose.
      ). However, the space between the OM and IM is restricted, and further study is required to determine how peptidoglycan and lignin interact.

      Conclusion

      Within a biorefinery context, where substrate concentrations are at the limit of what the organism can tolerate and metabolize, passive permeation for most LRCs through the OM will likely not be rate-limiting to bioconversion. Permeability coefficients for the OM are broadly similar to the IM counterparts, with the OM rate limiting for the most lipophilic LRCs while the IM is rate limiting for hydrophilic LRCs. The differences between the OM and the IM were largely attributable to the LPS leaflet found on the OM. Diffusion with the LPS layer is substantially slower than elsewhere in the OM, as interactions with LRCs and solution are made and broken with the LPS glycosylations. Charged substrates, however, permeate exceptionally slowly through the OM, limited by the unfavorable energetics for desolvating an ion as it passes through lipophilic membrane regions. Thus, we expect transport proteins to be intimately involved in carboxylate uptake by P. putida, while other LRCs permeate freely through a passive mechanism, depending on metabolic processes to generate the concentration gradient needed to drive flux across the OM.
      From a process design perspective, the ability for many LRCs to permeate passively eliminates some constraints on metabolic engineering approaches in other systems. For instance, sucrose synthesized photosynthetically depends on the overexpression of a sucrose–proton symporter for export at high titers (
      • Lin P.-C.
      • Zhang F.
      • Pakrasi H.B.
      Enhanced production of sucrose in the fast-growing cyanobacterium Synechococcus elongatus UTEX 2973.
      ), as sucrose cannot passively permeate bacterial membranes at an appreciable rate. While the final product from metabolically funneling LRCs may similarly require a specific transport protein for harvesting product, the inputs for bioconversion can be sourced without substantial engineering effort invested in importing the diverse chemistries represented by the range of LRCs tested. Thus, bacteria with flexible metabolisms that can process a wide range of LRCs and waste aromatics from other sources represent excellent candidates for future metabolic funneling without engineering specific transport proteins.

      Experimental procedures

      The general procedure for computing OM permeability closely follows prior permeation studies in the group for small molecules across biological membranes (
      • Vermaas J.V.
      • Dixon R.A.
      • Chen F.
      • Mansfield S.D.
      • Boerjan W.
      • Ralph J.
      • et al.
      Passive membrane transport of lignin-related compounds.
      ,
      • Vermaas J.V.
      • Beckham G.T.
      • Crowley M.F.
      Membrane permeability of fatty acyl compounds studied via molecular simulation.
      ,
      • Vermaas J.V.
      • Bentley G.J.
      • Beckham G.T.
      • Crowley M.F.
      Membrane permeability of terpenoids explored with molecular simulation.
      ). First, we carry out equilibrium MD simulations for a set of 42 LRCs (Fig. 1) in the context of a model bacterial OM from P. putida (Fig. 10), checking for spontaneous permeation as was observed for other small molecules through biological membranes (
      • Vermaas J.V.
      • Beckham G.T.
      • Crowley M.F.
      Membrane permeability of fatty acyl compounds studied via molecular simulation.
      ,
      • Vermaas J.V.
      • Bentley G.J.
      • Beckham G.T.
      • Crowley M.F.
      Membrane permeability of terpenoids explored with molecular simulation.
      ,
      • Odinokov A.
      • Ostroumov D.
      Structural degradation and swelling of lipid bilayer under the action of benzene.
      ). The equilibrium MD complements REUS simulations that allow for permeability to be calculated through the inhomogeneous solubility-diffusion model (
      • Marrink S.-J.
      • Berendsen H.J.C.
      Simulation of water transport through a lipid membrane.
      ,
      • Lee C.T.
      • Comer J.
      • Herndon C.
      • Leung N.
      • Pavlova A.
      • Swift R.V.
      • et al.
      Simulation-based approaches for determining membrane permeability of small compounds.
      ,
      • Vermaas J.V.
      • Beckham G.T.
      • Crowley M.F.
      Membrane permeability of fatty acyl compounds studied via molecular simulation.
      ,
      • Gupta R.
      • Sridhar D.B.
      • Rai B.
      Molecular dynamics simulation study of permeation of molecules through skin lipid bilayer.
      ,
      • Bennion B.J.
      • Be N.A.
      • McNerney M.W.
      • Lao V.
      • Carlson E.M.
      • Valdez C.A.
      • et al.
      Predicting a drug’s membrane permeability: a computational model validated with in vitro permeability assay data.
      ), even if no spontaneous permeation events are observed and the timescale for permeation is unknown. The inputs and selected outputs are available on Zenodo.
      The chosen LRCs for this study (Fig. 1) were drawn from a subset of the LRCs where permeation across a model P. putida IM was calculated (
      • Vermaas J.V.
      • Dixon R.A.
      • Chen F.
      • Mansfield S.D.
      • Boerjan W.
      • Ralph J.
      • et al.
      Passive membrane transport of lignin-related compounds.
      ), as the larger simulation systems required to capture a glycosylated OM made repeating the same suite of compounds cost prohibitive. The 42 selected LRCs in Figure 1 reflect the diverse chemistries seen in lignin but are tailored to focus on monomeric products that could result from native lignin polymer deconstruction (
      • Katahira R.
      • Mittal A.
      • McKinney K.
      • Chen X.
      • Tucker M.P.
      • Johnson D.K.
      • et al.
      Base-catalyzed depolymerization of biorefinery lignins.
      ,
      • Talebi Amiri M.
      • Dick G.R.
      • Questell-Santiago Y.M.
      • Luterbacher J.S.
      Fractionation of lignocellulosic biomass to produce uncondensed aldehyde-stabilized lignin.
      ). Other simple aromatics tested reflect species that might be encountered by a bacterial host in a mixed waste stream (
      • Borchert A.J.
      • Henson W.R.
      • Beckham G.T.
      Challenges and opportunities in biological funneling of heterogeneous and toxic substrates beyond lignin.
      ,
      • Black B.A.
      • Michener W.E.
      • Ramirez K.J.
      • Biddy M.J.
      • Knott B.C.
      • Jarvis M.W.
      • et al.
      Aqueous stream characterization from biomass fast pyrolysis and catalytic fast pyrolysis.
      ,
      • Henson W.R.
      • Meyers A.W.
      • Jayakody L.N.
      • DeCapite A.
      • Black B.A.
      • Michener W.E.
      • et al.
      Biological upgrading of pyrolysis-derived wastewater: engineering Pseudomonas putida for alkylphenol, furfural, and acetone catabolism and (methyl)muconic acid production.
      ).

      Membrane construction

      The OM of P. putida is composed of diverse LPSs). Lipidomics and glycomics have not elucidated the distribution of potential LPS species in enough detail for molecular modeling. Thus, our OM model uses 30 diverse LPS molecules modeled on what is found in Pseudomonas aeruginosa and makes specific structural choices to support our research objectives. The 30 LPS molecules for the OM were constructed through CHARMM-GUI (
      • Jo S.
      • Wu E.L.
      • Stuhlsatz D.
      • Klauda J.B.
      • MacKerell A.D.
      • Widmalm G.
      • et al.
      Lipopolysaccharide membrane building and simulation.
      ,
      • Wu E.L.
      • Cheng X.
      • Jo S.
      • Rui H.
      • Song K.C.
      • Dávila-Contreras E.M.
      • et al.
      CHARMM-GUI membrane builder toward realistic biological membrane simulations.
      ,
      • Jo S.
      • Kim T.
      • Iyer V.G.
      • Im W.
      CHARMM-GUI: a web-based graphical user interface for CHARMM.
      ,
      • Lee J.
      • Patel D.S.
      • Ståhle J.
      • Park S.-J.
      • Kern N.R.
      • Kim S.
      • et al.
      CHARMM-GUI membrane builder for complex biological membrane simulations with glycolipids and lipoglycans.
      ), six for each of the five available lipid A structures, using core 1B to represent the glycosylation pattern. Unlike other recent studies for the OM in P. aeruginosa (
      • Ĺopez C.A.
      • Zgurskaya H.
      • Gnanakaran S.
      Molecular characterization of the outer membrane of Pseudomonas aeruginosa.
      ), no O-antigen glycosylations were added to the LPS molecules. From a computational perspective, a minimal O-antigen would add at least another 4 nm of membrane thickness (
      • Ĺopez C.A.
      • Zgurskaya H.
      • Gnanakaran S.
      Molecular characterization of the outer membrane of Pseudomonas aeruginosa.
      ) for LRCs to traverse, and it is not immediately clear what the correct O-antigen length would be for P. putida. We estimate that additional glycosylation from a standard O-antigen would further increase the system size by around 40%. The increased system size decreases performance for individual simulations, as well as increasing the number of umbrellas that would be needed to span the bilayer when computing permeability. Thus, we estimate that adding in the O-antigen would effectively double the total computational cost for determining permeability. We approximate the impact for the O-antigen by assuming that the highest barrier we compute across the glycans is maintained for additional O-antigen thickness.
      The inner leaflet was also generated using CHARMM-GUI (
      • Wu E.L.
      • Cheng X.
      • Jo S.
      • Rui H.
      • Song K.C.
      • Dávila-Contreras E.M.
      • et al.
      CHARMM-GUI membrane builder toward realistic biological membrane simulations.
      ,
      • Jo S.
      • Kim T.
      • Iyer V.G.
      • Im W.
      CHARMM-GUI: a web-based graphical user interface for CHARMM.
      ,
      • Lee J.
      • Patel D.S.
      • Ståhle J.
      • Park S.-J.
      • Kern N.R.
      • Kim S.
      • et al.
      CHARMM-GUI membrane builder for complex biological membrane simulations with glycolipids and lipoglycans.
      ). The composition for the inner leaflet matches the composition for IM models for P. putida used previously (
      • Vermaas J.V.
      • Dixon R.A.
      • Chen F.
      • Mansfield S.D.
      • Boerjan W.
      • Ralph J.
      • et al.
      Passive membrane transport of lignin-related compounds.
      ), with an approximate 2:1 PE:PG headgroup ratio to match P. putida lipidomics studies (
      • Rühl J.
      • Hein E.-M.
      • Hayen H.
      • Schmid A.
      • Blank L.M.
      The glycerophospholipid inventory of Pseudomonas putida is conserved between strains and enables growth condition-related alterations.
      ). The final number of lipids in the inner leaflet (75) was chosen to match the membrane surface area for the LPS outer leaflet. The membrane is solvated and ionized within CHARMM-GUI, with 114 Ca2+, 114 K+, and 22 Cl ions. Calcium ions are chosen for ionization here based on the structural evidence for their importance to maintain OM stability (
      • Clifton L.A.
      • Skoda M.W.A.
      • Le Brun A.P.
      • Ciesielski F.
      • Kuzmenko I.
      • Holt S.A.
      • et al.
      Effect of divalent cation removal on the structure of gram-negative bacterial outer membrane models.
      ,
      • Clifton L.A.
      • Holt S.A.
      • Hughes A.V.
      • Daulton E.L.
      • Arunmanee W.
      • Heinrich F.
      • et al.
      An accurate in vitro model of the E. coli envelope.
      ). When the membrane geometry is complemented with small molecules, each molecular system is approximately 55,000 atoms in size, with equilibrated dimensions of approximately 75 × 75 × 90 Å.

      Equilibrium simulation

      To evaluate spontaneous permeability for each LRC considered (Fig. 1), 10 LRC molecules were added to the previously generated membrane (Fig. 10) to create independent simulation systems for each LRC. All 42 systems were simulated for 400 ns using NAMD 2.13 (
      • Phillips J.C.
      • Hardy D.J.
      • Maia J.D.C.
      • Stone J.E.
      • Ribeiro J.V.
      • Bernardi R.C.
      • et al.
      Scalable molecular dynamics on CPU and GPU architectures with NAMD.
      ). The systems were simulated using the CHARMM36 force field, including terms for lipids (
      • Klauda J.B.
      • Venable R.M.
      • Freites J.A.
      • O’Connor J.W.
      • Tobias D.J.
      • Mondragon-Ramirez C.
      • et al.
      Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types.
      ), carbohydrates (
      • Guvench O.
      • Greenr S.N.
      • Kamath G.
      • Brady J.W.
      • Venable R.M.
      • Pastor R.W.
      • et al.
      Additive empirical force field for hexopyranose monosaccharides.
      ,
      • Guvench O.
      • Hatcher E.
      • Venable R.M.
      • Pastor R.W.
      • MacKerell A.D.
      CHARMM additive all-atom force field for glycosidic linkages between hexopyranoses.
      ,
      • Guvench O.
      • Mallajosyula S.S.
      • Raman E.P.
      • Hatcher E.
      • Vanommeslaeghe K.
      • Foster T.J.
      • et al.
      CHARMM additive all-atom force field for carbohydrate derivatives and its utility in polysaccharide and carbohydrate-protein modeling.
      ), and lignin (
      • Vermaas J.V.
      • Petridis L.
      • Ralph J.
      • Crowley M.F.
      • Beckham G.T.
      Systematic parameterization of lignin for the CHARMM force field.
      ). The CHARMM36 General force field (
      • Vanommeslaeghe K.
      • Hatcher E.
      • Acharya C.
      • Kundu S.
      • Zhong S.
      • Shim J.
      • et al.
      CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields.
      ) was used for small organic molecules not covered by the lignin forcefield (Other Aromatics category in Fig. 1). Aqueous solution was represented with the TIP3 water model (
      • Jorgensen W.L.
      • Chandrasekhar J.
      • Madura J.D.
      • Impey R.W.
      • Klein M.L.
      • Pastor M.L.
      Comparison of simple potential functions for simulating liquid water.
      ) and ion parameters from Beglov and Roux (
      • Beglov D.
      • Roux B.
      Finite representation of an infinite bulk system: solvent boundary potential for computer simulations.
      ).
      Figure thumbnail gr10
      Figure 10Membrane construction snapshot. The heavy atoms for the lipid headgroups and acyl tails are shown as spheres, utilizing gray carbons to contrast with the yellow carbons featured in the glycosylations attached to the outer leaflet, represented with sticks. Other atoms have a consistent color scheme, with blue representing nitrogen, red representing oxygen, and brown representing phosphorus. Hydrogen and solvating water molecules, while present during simulation, are omitted for visual clarity.
      Other simulation parameters shared with the subsequent REUS calculations include a 12 Å nonbonded cutoff with switching applied past 10 Å and long-range electrostatics treatment using the particle mesh Ewald method (
      • Darden T.
      • York D.
      • Pedersen L.
      Particle mesh Ewald: an N log(N) method for Ewald sums in large systems.
      ,
      • Essmann U.
      • Perera L.
      • Berkowitz M.L.
      • Darden T.
      • Lee H.
      • Pedersen L.G.
      A smooth particle mesh Ewald method.
      ) with 1.2 Å grid spacing. Semi-isotropic pressure coupling was maintained via the Langevin piston method (
      • Feller S.E.
      • Zhang Y.
      • Pastor R.W.
      • Brooks B.R.
      Constant pressure molecular dynamics simulation: the Langevin piston method.
      ,
      • Martyna G.J.
      • Tobias D.J.
      • Klein M.L.
      Constant pressure molecular dynamics algorithms.
      ) to a pressure of 1 atm, with periodic cell dimensions along the membrane normal axis decoupled from growth along the membrane parallel axes. A Langevin thermostat (
      • Kubo R.
      The fluctuation-dissipation theorem.
      ) using γ=1 ps−1 maintained the system temperature at 310K. Each simulation timestep was 2 fs, enabled by using the SETTLE algorithm (
      • Miyamoto S.
      • Kollman P.a.
      Settle: an analytical version of the SHAKE and RATTLE algorithm for rigid water models.
      ) to fix bond lengths to hydrogen.

      Replica exchange umbrella sampling

      To obtain a free energy profile and local diffusivity for our selected LRCs along the reaction coordinate representing passive membrane transit across the OM, we carried out a set of REUS (
      • Sugita Y.
      • Kitao A.
      • Okamoto Y.
      Multidimensional replica-exchange method for free-energy calculations.
      ) (also called Hamiltonian replica exchange (
      • Fukunishi H.
      • Watanabe O.
      • Takada S.
      On the Hamiltonian replica exchange method for efficient sampling of biomolecular systems: application to protein structure prediction.
      ), window exchange umbrella sampling (
      • Park S.
      • Kim T.
      • Im W.
      Transmembrane helix assembly by window exchange umbrella sampling.
      ), or bias exchange umbrella sampling (
      • Moradi M.
      • Tajkhorshid E.
      Computational recipe for efficient description of large-scale conformational changes in biomolecular systems.
      )) simulations. Similar to conventional umbrella sampling, REUS samples high energy regions of the reaction coordinate range by applying a harmonic bias (umbrella) to the system, forcing uniform sampling across a reaction coordinate (
      • Torrie G.
      • Valleau J.
      Nonphysical sampling distributions in Monte Carlo free-energy estimation: umbrella sampling.
      ), rather than just the low energy regions of the free energy surface as would occur without an applied bias. Given known biases, unbiased free energy surfaces can be estimated from these biased simulations using well-established techniques (
      • Kumar S.
      • Rosenberg J.M.
      • Bouzida D.
      • Swendsen R.H.
      • Kollman P.A.
      The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method.
      ,
      • Ferguson A.L.
      BayesWHAM: a bayesian approach for free energy estimation, reweighting, and uncertainty quantification in the weighted histogram analysis method.
      ).
      For our system, the chosen reaction coordinate reflects permeation for the center of mass for a single LRC molecule across the membrane, going from -35Å on the inner leaflet side to 65Å on the outer leaflet side. As each umbrella restricts the sampling to a small fraction of the total reaction coordinate, a series of 180 independent umbrellas with 4 kcalmol−1 force constants were used to cover the entire reaction coordinate, yielding a 0.55 Å umbrella spacing. This spacing and force constant are consistent with what has worked previously to evaluate permeation in other systems (
      • Vermaas J.V.
      • Dixon R.A.
      • Chen F.
      • Mansfield S.D.
      • Boerjan W.
      • Ralph J.
      • et al.
      Passive membrane transport of lignin-related compounds.
      ,
      • Lee C.T.
      • Comer J.
      • Herndon C.
      • Leung N.
      • Pavlova A.
      • Swift R.V.
      • et al.
      Simulation-based approaches for determining membrane permeability of small compounds.
      ).
      The primary advantage of REUS relative to conventional umbrella sampling is ensuring contiguous sampling, as the individual umbrella biases are exchanged periodically according to a Metropolis criterion (
      • Fukunishi H.
      • Watanabe O.
      • Takada S.
      On the Hamiltonian replica exchange method for efficient sampling of biomolecular systems: application to protein structure prediction.
      ,
      • Sugita Y.
      • Kitao A.
      • Okamoto Y.
      Multidimensional replica-exchange method for free-energy calculations.
      ,
      • Metropolis N.
      • Rosenbluth A.W.
      • Rosenbluth M.N.
      • Teller A.H.
      • Teller E.
      Equation of state calculations by fast computing machines.
      ), accelerating sampling along the chosen reaction coordinate. Adjacent windows were attempted to be exchanged with alternating neighbors every 1 ps, consistent with advice to exchange often during replica exchange simulations (
      • Sindhikara D.J.
      • Emerson D.J.
      • Roitberg A.E.
      Exchange often and properly in replica exchange molecular dynamics.
      ). To set up a diverse pool of starting configurations for each replica, initial configurations for single LRC molecules in an OM system are drawn at random from the computed equilibrium trajectories. For computing permeabilities, REUS achieves the same level of accuracy as conventional umbrella sampling with shorter simulation times (
      • Lee C.T.
      • Comer J.
      • Herndon C.
      • Leung N.
      • Pavlova A.
      • Swift R.V.
      • et al.
      Simulation-based approaches for determining membrane permeability of small compounds.
      ), and based on past experience (
      • Vermaas J.V.
      • Dixon R.A.
      • Chen F.
      • Mansfield S.D.
      • Boerjan W.
      • Ralph J.
      • et al.
      Passive membrane transport of lignin-related compounds.
      ), each window was sampled for 40 ns.

      ABF calculation

      To assess the impact for the O-antigen on LRC permeation, we determined the local free energy and diffusivity profiles using ABF (
      • Darve E.
      • Pohorille A.
      Calculating free energies using average force.
      ,
      • Comer J.
      • Gumbart J.C.
      • Hénin J.
      • Lelièvre T.
      • Pohorille A.
      • Chipot C.
      The adaptive biasing force method: everything you always wanted to know but were afraid to ask.
      ) simulations. The membrane setup matches the membrane composition described previously but adds O-antigen consistent with the P. aeruginosa default composition from the CHARMM-GUI membrane builder (
      • Lee J.
      • Patel D.S.
      • Ståhle J.
      • Park S.-J.
      • Kern N.R.
      • Kim S.
      • et al.
      CHARMM-GUI membrane builder for complex biological membrane simulations with glycolipids and lipoglycans.
      ). This adds 15 rhamnose monomers to the membrane, increasing the height for the membrane considerably (Fig. 11). Once solvated and ionized as noted above, the initial system is 60 × 60 × 150 Å.
      Figure thumbnail gr11
      Figure 11Initial snapshot for the membrane system with the O-antigen added. The representations used are analogous to , with the 15 additional rhamnose monomers on each LPS head shown with green carbons, while the glycosylated core retains the yellow carbons from . LPS, lipopolysaccharide.
      To sample the membrane normal dimension efficiently, the full membrane normal span was split into 15 independently sampled reaction coordinates. Each collective variable defining the space was 13 Å wide with 2 Å overlap between adjacent collective variables, spanning from -40 Å below the membrane midplane to 127 Å above the membrane midplane. A single guaiacol molecule sampled each bounded reaction coordinate, such that 15 guaiacol molecules in total were present within the membrane-water system. To further increase sampling, eight ABF simulations were performed simultaneously and shared biasing information through the multiple walker framework every 5000 steps (
      • Minoukadeh K.
      • Chipot C.
      • Lelièvre T.
      Potential of mean force calculations: a multiple-walker adaptive biasing force approach.
      ) within NAMD 2.14 (
      • Phillips J.C.
      • Hardy D.J.
      • Maia J.D.C.
      • Stone J.E.
      • Ribeiro J.V.
      • Bernardi R.C.
      • et al.
      Scalable molecular dynamics on CPU and GPU architectures with NAMD.
      ). Thus, the eight individual 378 ns simulations performed represent over 45μs of sampling for the membrane normal reaction coordinate, substantially longer than the 7.2μs needed for a converged PMF with REUS.

      Analysis

      The analysis procedures largely follow our prior work on permeability across biological membranes (
      • Vermaas J.V.
      • Dixon R.A.
      • Chen F.
      • Mansfield S.D.
      • Boerjan W.
      • Ralph J.
      • et al.
      Passive membrane transport of lignin-related compounds.
      ,
      • Vermaas J.V.
      • Beckham G.T.
      • Crowley M.F.
      Membrane permeability of fatty acyl compounds studied via molecular simulation.
      ,
      • Vermaas J.V.
      • Bentley G.J.
      • Beckham G.T.
      • Crowley M.F.
      Membrane permeability of terpenoids explored with molecular simulation.
      ). Python scripts combining VMD (
      • Humphrey W.
      • Dalke A.
      • Schulten K.
      VMD: visual molecular dynamics.
      ), numpy (
      • van der Walt S.
      • Colbert S.C.
      • Varoquaux G.
      The NumPy array: a structure for efficient numerical computation.
      ,
      • Harris C.R.
      • Millman K.J.
      • van der Walt S.J.
      • Gommers R.
      • Virtanen P.
      • Cournapeau D.
      • et al.
      Array programming with NumPy.
      ), and scipy quantified simulation trajectory observables. Molecular visualizations were rendered using VMD (
      • Humphrey W.
      • Dalke A.
      • Schulten K.
      VMD: visual molecular dynamics.
      ) with its built-in GPU-accelerated OptiX raytracer (
      • Stone J.E.
      • Vandivort K.L.
      • Schulten K.
      GPU-accelerated molecular visualization on petascale supercomputing platforms.
      ), and dataplots were made with matplotlib (
      • Hunter J.D.
      Matplotlib: a 2D graphics environment.
      ). From the equilibrium trajectories, we track the location of individual LRC molecules to identify complete crossing events, as well as membrane-depth dependent estimates for solvent diffusion. Diffusion for water and lipids within the equilibrium trajectories are determined via the Einstein relation (
      • Einstein A.
      Uber die von der molekularkinetischen Theorie der Wärmë geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen.
      ) in either two dimensions for lateral diffusion or in one dimension for diffusion along the membrane normal axis. We do so by tracking displacements at 20 ps intervals and binning the resulting displacements based on position to estimate position-dependent water diffusion from equilibrium simulation. The carbon-deuterium order parameter measuring structure for the acyl tails (-SCD) through NMR is approximated by the carbon-hydrogen order parameter used in parameterizing the lipid force field (
      • Klauda J.B.
      • Venable R.M.
      • Freites J.A.
      • O’Connor J.W.
      • Tobias D.J.
      • Mondragon-Ramirez C.
      • et al.
      Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types.
      ). The REUS calculations are the basis for calculating the permeability and partitioning coefficients, based on free energy and diffusivity profiles derived from the REUS trajectories.
      The free energies G and diffusivities D along a reaction coordinate ξ bounded by upper and lower boundaries ξu and ξl are converted into permeability coefficients (Pm) through the inhomogeneous solubility-diffusion model (
      • Marrink S.-J.
      • Berendsen H.J.C.
      Simulation of water transport through a lipid membrane.
      ,
      • Lee C.T.
      • Comer J.
      • Herndon C.
      • Leung N.
      • Pavlova A.
      • Swift R.V.
      • et al.
      Simulation-based approaches for determining membrane permeability of small compounds.
      ) originally postulated by Diamond and Katz (
      • Diamond J.M.
      • Katz Y.
      Interpretation of nonelectrolyte partition coefficients between dimyristoyl lecithin and water.
      ).
      Pm=[ξlξuexp(ΔG(ξ)β)D(ξ)dξ]1
      (3)


      We further break down the permeability coefficient into three subcomponents representing a compound crossing the hydrophobic region (Pmc), exiting the membrane across the glycosylations (Pmg), and exiting the membrane across the unglycosylated side (Pmu) by splitting the integration range within Equation 3 into three parts. The integration boundaries splitting the integral to create the three subcomponents were determined by the free energy minima in the range (-22,-5) for the unglycosylated side boundary and the range (5, 22) on the glycosylated side boundary. The partition coefficient depends solely on the free energy difference between aqueous solution and the membrane.
      logP=GaqGmembraneRTln10=ΔGpartitionRTln10
      (4)


      In this form, we use the simple approximation that Gmembrane is the lower free energy minimum at the integration boundaries determined for permeability. The free energy reference state while integrating Equation 3 to obtain the subcomponents is chosen to be the minimum free energy. As a result, the permeability coefficient for traversing the entire membrane depends on the partition coefficient computed from Equation 4 to reset the free energy reference to capture the real kinetics for moving from solution to solution. In mathematical terms, this works out to be:
      log[Pm]=logP+log[(Pmu1+Pmc1+Pmg1)1]
      (5)


      To estimate the O-antigen thickness needed to decrease permeability by tenfold, we first find the position ξmax where the free energy is maximized within the LPS layer. At this position, the integrand for Equation 3 (exp(ΔG(ξmax)β)D(ξmax)) is maximized, which represents the incremental resistance to permeation (
      • Marrink S.-J.
      • Berendsen H.J.C.
      Simulation of water transport through a lipid membrane.
      ,
      • Lee C.T.
      • Comer J.
      • Herndon C.
      • Leung N.
      • Pavlova A.
      • Swift R.V.
      • et al.
      Simulation-based approaches for determining membrane permeability of small compounds.
      ) as the membrane thickens due to longer/larger O-antigen bands. We solve for the thickness L needed to increase the resistivity R = Pm−1 ten-fold, which would reduce the permeability by an order of magnitude.
      10R=exp(ΔG(ξmax)β)D(ξmax)L
      (6)


      The free energy profile from the REUS trajectories were calculated from a modified version of BayesWHAM (
      • Ferguson A.L.
      BayesWHAM: a bayesian approach for free energy estimation, reweighting, and uncertainty quantification in the weighted histogram analysis method.
      ), which uses Gibbs sampling of the known Dirichlet prior (
      • Habeck M.
      Bayesian estimation of free energies from equilibrium simulations.
      ) to rapidly assess uncertainty. We use the free energy in solution as our zero point and force the equivalence for the free energy at the edges of our reaction coordinate by treating the reaction coordinate as periodic. The periodicity assumption has the largest impact on the calculated profiles for charged LRCs, namely the benzoates and cinnamates, as dehydrating ionic species during permeation can yield large sampling errors that are remedied in part by symmetrization (
      • Sun L.
      • Bertelshofer F.
      • Greiner G.
      • Böckmann R.A.
      Characteristics of sucrose transport through the sucrose-specific porin ScrY studied by molecular dynamics simulations.
      ,
      • Golla V.K.
      • Prajapati J.D.
      • Joshi M.
      • Kleinekathöfer U.
      Exploration of free energy surfaces across a membrane channel using metadynamics and umbrella sampling.
      ). Diffusivity profiles for Equation 3 were computed directly from the variance and autocorrelation time of the biased motion along the reaction coordinate (
      • Hummer G.
      Position-dependent diffusion coefficients and free energies from Bayesian analysis of equilibrium and replica molecular dynamics simulations.
      ), which is less sensitive to the harmonic restraint force than alternative calculation approaches (
      • Gaalswyk K.
      • Awoonor-Williams E.
      • Rowley C.N.
      Generalized Langevin methods for calculating transmembrane diffusivity.
      ) and is not influenced by concerns surrounding momentum removal (
      • Fujimoto K.
      • Nagai T.
      • Yamaguchi T.
      Momentum removal to obtain the position-dependent diffusion constant in constrained molecular dynamics simulation.
      ). To compute auto-correlation times needed for the diffusivity calculation, we fit the autocorrelation function to an exponentially decaying function within the 1 ps intervals between exchanges. This has been previously identified to be a suitable method for estimating small molecule diffusivity from REUS simulations (
      • Vermaas J.V.
      • Beckham G.T.
      • Crowley M.F.
      Membrane permeability of fatty acyl compounds studied via molecular simulation.
      ). The free energy profile from the ABF simulation was stitched together from the 15 independent reaction coordinates in numpy (
      • van der Walt S.
      • Colbert S.C.
      • Varoquaux G.
      The NumPy array: a structure for efficient numerical computation.
      ,
      • Harris C.R.
      • Millman K.J.
      • van der Walt S.J.
      • Gommers R.
      • Virtanen P.
      • Cournapeau D.
      • et al.
      Array programming with NumPy.
      ). Diffusivity estimates directly from the ABF simulation were computed using DiffusionFusion, a C implementation for calculating position-dependent diffusivity from ABF calculations (
      • Comer J.
      • Chipot C.
      • González-Nilo F.D.
      Calculating position-dependent diffusivity in biased molecular dynamics simulations.
      ).

      Data availability

      The reduced directory structure that includes analysis scripts, inputs, and selected raw outputs used for this publication is available from http://doi.org/10.5281/zenodo.5794252.
      The complete directory structure is available upon request.

      Supporting information

      The Supporting Information contains a pdf file with Table S1 reporting the partitioning and permeability coefficients for all 42 compounds. Figures S1–S41 supplement Figure 2. Figures S42–S51 report the free energy and diffusivity profiles for all 42 compounds, analogous to Figure 5. The free energy and diffusivity profiles are also tabulated numerically as an excel spreadsheet. Animation S1 shows the permeation of a single syringol molecule across the membrane in unbiased simulation, complementing Figure 2.

      Conflict of interest

      The authors declare that they have no conflicts of interest with the contents of this article.

      Acknowledgments

      We would like to thank an anonymous reviewer from our previous work (
      • Vermaas J.V.
      • Dixon R.A.
      • Chen F.
      • Mansfield S.D.
      • Boerjan W.
      • Ralph J.
      • et al.
      Passive membrane transport of lignin-related compounds.
      ) for wondering how LRC permeation would differ between the IM and OM. This work was authored in part by Alliance for Sustainable Energy, LLC, the manager and operator of the National Renewable Energy Laboratory for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. This research used resources of the 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 project made use of the high performance computing resources provided by the NREL Computational Sciences Center, which is supported by the DOE Office of EERE under Contract No. DE-AC36-08GO28308 . This work was supported in part through computational resources and services provided by the Institute for Cyber-Enabled Research at Michigan State University .

      Author contributions

      J. V. V., M. F. C., and G. T. B. conceptualization; J. V. V. investigation; J. V. V. writing–original draft; M. F. C. and G. T. B. writing–review and editing.

      Funding and additional information

      J. V. V., M. F. C., and G. T. B. acknowledge funding from the U.S. Department of Energy Office of Energy Efficiency and Renewable Energy Bioenergy Technologies Office. G. T. B. also acknowledges funding from the Center for Bioenergy Innovation, which is a U.S. DOE Bioenergy Research Center supported by the Office of Biological and Environmental Research in the DOE Office of Science.

      Supporting information

      References

        • Boerjan W.
        • Ralph J.
        • Baucher M.
        Lignin biosynthesis.
        Annu. Rev. Plant Biol. 2003; 54: 519-546
        • Ragauskas A.J.
        • Beckham G.T.
        • Biddy M.J.
        • Chandra R.
        • Chen F.
        • Davis M.F.
        • et al.
        Lignin valorization: improving lignin processing in the biorefinery.
        Science. 2014; 344: 1246843
        • Janusz G.
        • Pawlik A.
        • Sulej J.
        • Świderska-Burek U.
        • Jarosz-Wilkołazka A.
        • Paszczyński A.
        Lignin degradation: microorganisms, enzymes involved, genomes analysis and evolution.
        FEMS Microbiol. Rev. 2017; 41: 941-962
        • Bugg T.D.H.
        • Ahmad M.
        • Hardiman E.M.
        • Rahmanpour R.
        Pathways for degradation of lignin in bacteria and fungi.
        Nat. Prod. Rep. 2011; 28: 1883
        • Borchert A.J.
        • Henson W.R.
        • Beckham G.T.
        Challenges and opportunities in biological funneling of heterogeneous and toxic substrates beyond lignin.
        Curr. Opin. Biotechnol. 2022; 73: 1-13
        • Linger J.G.
        • Vardon D.R.
        • Guarnieri M.T.
        • Karp E.M.
        • Hunsinger G.B.
        • Franden M.A.
        • et al.
        Lignin valorization through integrated biological funneling and chemical catalysis.
        Proc. Natl. Acad. Sci. U. S. A. 2014; 111: 12013-12018
        • Becker J.
        • Wittmann C.
        From systems biology to metabolically engineered cells - an omics perspective on the development of industrial microbes.
        Curr. Opin. Microbiol. 2018; 45: 180-188
        • Abdelaziz O.Y.
        • Brink D.P.
        • Prothmann J.
        • Ravi K.
        • Sun M.
        • Garćıa-Hidalgo J.
        • et al.
        Biological valorization of low molecular weight lignin.
        Biotechnol. Adv. 2016; 34: 1318-1346
        • Kamimura N.
        • Takahashi K.
        • Mori K.
        • Araki T.
        • Fujita M.
        • Higuchi Y.
        • et al.
        Bacterial catabolism of lignin-derived aromatics: new findings in a recent decade: update on bacterial lignin catabolism.
        Environ. Microbiol. Rep. 2017; 9: 679-705
        • Erickson E.
        • Bleem A.
        • Kuatsjah E.
        • Werner A.Z.
        • DuBois J.L.
        • McGeehan J.E.
        • et al.
        Critical enzyme reactions in aromatic catabolism for microbial lignin conversion.
        Nat. Catal. 2022; 5: 86-98
        • Beckham G.T.
        • Johnson C.W.
        • Karp E.M.
        • Salvachúa D.
        • Vardon D.R.
        Opportunities and challenges in biological lignin valorization.
        Curr. Opin. Biotechnol. 2016; 42: 40-53
        • Johnson C.W.
        • Salvachúa D.
        • Rorrer N.A.
        • Black B.A.
        • Vardon D.R.
        • St. John P.C.
        • et al.
        Innovative chemicals and materials from bacterial aromatic catabolic pathways.
        Joule. 2019; 3: 1523-1537
        • del Cerro C.
        • Erickson E.
        • Dong T.
        • Wong A.R.
        • Eder E.K.
        • Purvine S.O.
        • et al.
        Intracellular pathways for lignin catabolism in white-rot fungi.
        Proc. Natl. Acad. Sci. U. S. A. 2021; 118e2017381118
        • Silhavy T.J.
        • Kahne D.
        • Walker S.
        The bacterial cell envelope.
        Cold Spring Harb. Perspect. Biol. 2010; 2a000414
        • Vermaas J.V.
        • Dixon R.A.
        • Chen F.
        • Mansfield S.D.
        • Boerjan W.
        • Ralph J.
        • et al.
        Passive membrane transport of lignin-related compounds.
        Proc. Natl. Acad. Sci. U. S. A. 2019; 116: 23117-23123
        • Salvachúa D.
        • Werner A.Z.
        • Pardo I.
        • Michalska M.
        • Black B.A.
        • Donohoe B.S.
        • et al.
        Outer membrane vesicles catabolize lignin-derived aromatic compounds in Pseudomonas putida KT2440.
        Proc. Natl. Acad. Sci. U. S. A. 2020; 117: 9302-9310
        • Perkins M.L.
        • Schuetz M.
        • Unda F.
        • Chen K.T.
        • Bally M.B.
        • Kulkarni J.A.
        • et al.
        Monolignol export by diffusion down a polymerization-induced concentration gradient.
        Plant Cell. 2022; 34: 2080-2095
        • Martinotti C.
        • Ruiz-Perez L.
        • Deplazes E.
        • Mancera R.L.
        Molecular dynamics simulation of small molecules interacting with biological membranes.
        ChemPhysChem. 2020; 21: 1486-1514
        • Venable R.M.
        • Krämer A.
        • Pastor R.W.
        Molecular dynamics simulations of membrane permeability.
        Chem. Rev. 2019; 119: 5954-5997
        • Jo S.
        • Wu E.L.
        • Stuhlsatz D.
        • Klauda J.B.
        • MacKerell A.D.
        • Widmalm G.
        • et al.
        Lipopolysaccharide membrane building and simulation.
        Met. Mol. Biol. 2015; 1273: 391-406
        • Khalid S.
        • Piggot T.J.
        • Samsudin F.
        Atomistic and coarse grain simulations of the cell envelope of gram-negative bacteria: what have we learned?.
        Acc. Chem. Res. 2019; 52: 180-188
        • Ĺopez C.A.
        • Zgurskaya H.
        • Gnanakaran S.
        Molecular characterization of the outer membrane of Pseudomonas aeruginosa.
        Biochim. Biophys. Acta - Biomembr. 2020; 1862183151
        • Piggot T.J.
        • Holdbrook D.A.
        • Khalid S.
        Electroporation of the E. coli and S. Aureus membranes: molecular dynamics simulations of complex bacterial membranes.
        J. Phys. Chem. B. 2011; 115: 13381-13388
        • Soares T.
        • Straatsma T.
        Assessment of the convergence of molecular dynamics simulations of lipopolysaccharide membranes.
        Mol. Simul. 2008; 34: 295-307
        • Balusek C.
        • Gumbart J.C.
        Role of the native outer-membrane environment on the transporter BtuB.
        Biophys. J. 2016; 111: 1409-1417
        • Luna E.
        • Kim S.
        • Gao Y.
        • Widmalm G.
        • Im W.
        Influences of Vibrio cholerae lipid A types on LPS bilayer properties.
        J. Phys. Chem. B. 2021; 125: 2105-2112
        • Wu E.L.
        • Cheng X.
        • Jo S.
        • Rui H.
        • Song K.C.
        • Dávila-Contreras E.M.
        • et al.
        CHARMM-GUI membrane builder toward realistic biological membrane simulations.
        J. Comput. Chem. 2014; 35: 1997-2004
        • Piggot T.J.
        • Allison J.R.
        • Sessions R.B.
        • Essex J.W.
        On the calculation of acyl chain order parameters from lipid simulations.
        J. Chem. Theor. Comput. 2017; 13: 5683-5696
        • Marrink S.-J.
        • Berendsen H.J.C.
        Simulation of water transport through a lipid membrane.
        J. Phys. Chem. 1994; 98: 4155-4168
        • Lee C.T.
        • Comer J.
        • Herndon C.
        • Leung N.
        • Pavlova A.
        • Swift R.V.
        • et al.
        Simulation-based approaches for determining membrane permeability of small compounds.
        J. Chem. Inf. Model. 2016; 56: 721-733
        • Ghysels A.
        • Venable R.M.
        • Pastor R.W.
        • Hummer G.
        Position-dependent diffusion tensors in anisotropic media from simulation: oxygen transport in and through membranes.
        J. Chem. Theor. Comput. 2017; 13: 2962-2976
        • De Vos O.
        • Venable R.M.
        • Van Hecke T.
        • Hummer G.
        • Pastor R.W.
        • Ghysels A.
        Membrane permeability: characteristic times and lengths for oxygen and a simulation-based test of the inhomogeneous solubility-diffusion model.
        J. Chem. Theor. Comput. 2018; 14: 3811-3824
        • Sicard F.
        • Koskin V.
        • Annibale A.
        • Rosta E.
        Position-dependent diffusion from biased simulations and markov state model analysis.
        J. Chem. Theor. Comput. 2021; 17: 2022-2033
        • Wilson A.N.
        • Dutta A.
        • Black B.A.
        • Mukarakate C.
        • Magrini K.
        • Schaidle J.A.
        • et al.
        Valorization of aqueous waste streams from thermochemical biorefineries.
        Green. Chem. 2019; 21: 4217-4230
        • Wynands B.
        • Lenzen C.
        • Otto M.
        • Koch F.
        • Blank L.M.
        • Wierckx N.
        Metabolic engineering of Pseudomonas taiwanensis VLB120 with minimal genomic modifications for high-yield phenol production.
        Metab. Eng. 2018; 47: 121-133
        • Wierckx N.J.P.
        • Ballerstedt H.
        • de Bont J.A.M.
        • Wery J.
        Engineering of solvent-tolerant Pseudomonas putida S12 for bioproduction of phenol from glucose.
        Appl. Environ. Microbiol. 2005; 71: 8221-8227
        • Mallinson S.J.B.
        • Machovina M.M.
        • Silveira R.L.
        • Garcia-Borra`s M.
        • Gallup N.
        • Johnson C.W.
        • et al.
        A promiscuous cytochrome P450 aromatic O-demethylase for lignin bioconversion.
        Nat. Commun. 2018; 9: 2487
        • Notonier S.
        • Werner A.Z.
        • Kuatsjah E.
        • Dumalo L.
        • Abraham P.E.
        • Hatmaker E.A.
        • et al.
        Metabolism of syringyl lignin-derived compounds in Pseudomonas putida enables convergent production of 2-pyrone-4,6-dicarboxylic acid.
        Metab. Eng. 2021; 65: 111-122
        • de Vries L.
        • MacKay H.A.
        • Smith R.A.
        • Mottiar Y.
        • Karlen S.D.
        • Unda F.
        • et al.
        HBMT1, a BAHD-family monolignol acyltransferase, mediates lignin acylation in poplar.
        Plant Physiol. 2022; 188: 1014-1027https://doi.org/10.1093/plphys/kiab546
        • Carpenter T.S.
        • Parkin J.
        • Khalid S.
        The free energy of small solute permeation through the Escherichia coli outer membrane has a distinctly asymmetric profile.
        J. Phys. Chem. Lett. 2016; 7: 3446-3451
        • Gumbart J.C.
        • Beeby M.
        • Jensen G.J.
        • Roux B.
        Escherichia coli peptidoglycan structure and mechanics as predicted by atomic-scale simulations.
        PLoS Comput. Biol. 2014; 10e1003475
        • Vermaas J.V.
        • Petridis L.
        • Qi X.
        • Schulz R.
        • Lindner B.
        • Smith J.C.
        Mechanism of lignin inhibition of enzymatic biomass deconstruction.
        Biotechnol. Biofuels. 2015; 8: 217
        • Vermaas J.V.
        • Crowley M.F.
        • Beckham G.T.
        A quantitative molecular atlas for interactions between lignin and cellulose.
        ACS Sustain. Chem. Eng. 2019; 7: 19570-19583
        • Lin P.-C.
        • Zhang F.
        • Pakrasi H.B.
        Enhanced production of sucrose in the fast-growing cyanobacterium Synechococcus elongatus UTEX 2973.
        Sci. Rep. 2020; 10: 390
        • Vermaas J.V.
        • Beckham G.T.
        • Crowley M.F.
        Membrane permeability of fatty acyl compounds studied via molecular simulation.
        J. Phys. Chem. B. 2017; 121: 11311-11324
        • Vermaas J.V.
        • Bentley G.J.
        • Beckham G.T.
        • Crowley M.F.
        Membrane permeability of terpenoids explored with molecular simulation.
        J. Phys. Chem. B. 2018; 122: 10349-10361
        • Odinokov A.
        • Ostroumov D.
        Structural degradation and swelling of lipid bilayer under the action of benzene.
        J. Phys. Chem. B. 2015; 119: 15006-15013
        • Gupta R.
        • Sridhar D.B.
        • Rai B.
        Molecular dynamics simulation study of permeation of molecules through skin lipid bilayer.
        J. Phys. Chem. B. 2016; 120: 8987-8996
        • Bennion B.J.
        • Be N.A.
        • McNerney M.W.
        • Lao V.
        • Carlson E.M.
        • Valdez C.A.
        • et al.
        Predicting a drug’s membrane permeability: a computational model validated with in vitro permeability assay data.
        J. Phys. Chem. B. 2017; 121: 5228-5237
        • Katahira R.
        • Mittal A.
        • McKinney K.
        • Chen X.
        • Tucker M.P.
        • Johnson D.K.
        • et al.
        Base-catalyzed depolymerization of biorefinery lignins.
        ACS Sustain. Chem. Eng. 2016; 4: 1474-1486
        • Talebi Amiri M.
        • Dick G.R.
        • Questell-Santiago Y.M.
        • Luterbacher J.S.
        Fractionation of lignocellulosic biomass to produce uncondensed aldehyde-stabilized lignin.
        Nat. Protoc. 2019; 14: 921-954
        • Black B.A.
        • Michener W.E.
        • Ramirez K.J.
        • Biddy M.J.
        • Knott B.C.
        • Jarvis M.W.
        • et al.
        Aqueous stream characterization from biomass fast pyrolysis and catalytic fast pyrolysis.
        ACS Sustain. Chem. Eng. 2016; 4: 6815-6827
        • Henson W.R.
        • Meyers A.W.
        • Jayakody L.N.
        • DeCapite A.
        • Black B.A.
        • Michener W.E.
        • et al.
        Biological upgrading of pyrolysis-derived wastewater: engineering Pseudomonas putida for alkylphenol, furfural, and acetone catabolism and (methyl)muconic acid production.
        Metab. Eng. 2021; 68: 14-25
        • Jo S.
        • Kim T.
        • Iyer V.G.
        • Im W.
        CHARMM-GUI: a web-based graphical user interface for CHARMM.
        J. Comput. Chem. 2008; 29: 1859-1865
        • Lee J.
        • Patel D.S.
        • Ståhle J.
        • Park S.-J.
        • Kern N.R.
        • Kim S.
        • et al.
        CHARMM-GUI membrane builder for complex biological membrane simulations with glycolipids and lipoglycans.
        J. Chem. Theor. Comput. 2019; 15: 775-786
        • Rühl J.
        • Hein E.-M.
        • Hayen H.
        • Schmid A.
        • Blank L.M.
        The glycerophospholipid inventory of Pseudomonas putida is conserved between strains and enables growth condition-related alterations.
        Microb. Biotechnol. 2012; 5: 45-58
        • Clifton L.A.
        • Skoda M.W.A.
        • Le Brun A.P.
        • Ciesielski F.
        • Kuzmenko I.
        • Holt S.A.
        • et al.
        Effect of divalent cation removal on the structure of gram-negative bacterial outer membrane models.
        Langmuir. 2015; 31: 404-412
        • Clifton L.A.
        • Holt S.A.
        • Hughes A.V.
        • Daulton E.L.
        • Arunmanee W.
        • Heinrich F.
        • et al.
        An accurate in vitro model of the E. coli envelope.
        Angew. Chem. Int. Ed. 2015; 54: 11952-11955
        • Phillips J.C.
        • Hardy D.J.
        • Maia J.D.C.
        • Stone J.E.
        • Ribeiro J.V.
        • Bernardi R.C.
        • et al.
        Scalable molecular dynamics on CPU and GPU architectures with NAMD.
        J. Chem. Phys. 2020; 153044130
        • Klauda J.B.
        • Venable R.M.
        • Freites J.A.
        • O’Connor J.W.
        • Tobias D.J.
        • Mondragon-Ramirez C.
        • et al.
        Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types.
        J. Phys. Chem. B. 2010; 114: 7830-7843
        • Guvench O.
        • Greenr S.N.
        • Kamath G.
        • Brady J.W.
        • Venable R.M.
        • Pastor R.W.
        • et al.
        Additive empirical force field for hexopyranose monosaccharides.
        J. Comput. Chem. 2008; 29: 2543-2564
        • Guvench O.
        • Hatcher E.
        • Venable R.M.
        • Pastor R.W.
        • MacKerell A.D.
        CHARMM additive all-atom force field for glycosidic linkages between hexopyranoses.
        J. Chem. Theor. Comput. 2009; 5: 2353-2370
        • Guvench O.
        • Mallajosyula S.S.
        • Raman E.P.
        • Hatcher E.
        • Vanommeslaeghe K.
        • Foster T.J.
        • et al.
        CHARMM additive all-atom force field for carbohydrate derivatives and its utility in polysaccharide and carbohydrate-protein modeling.
        J. Chem. Theor. Comput. 2011; 7: 3162-3180
        • Vermaas J.V.
        • Petridis L.
        • Ralph J.
        • Crowley M.F.
        • Beckham G.T.
        Systematic parameterization of lignin for the CHARMM force field.
        Green. Chem. 2019; 21: 109-122
        • Vanommeslaeghe K.
        • Hatcher E.
        • Acharya C.
        • Kundu S.
        • Zhong S.
        • Shim J.
        • et al.
        CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields.
        J. Comput. Chem. 2010; 31: 671-690
        • Jorgensen W.L.
        • Chandrasekhar J.
        • Madura J.D.
        • Impey R.W.
        • Klein M.L.
        • Pastor M.L.
        Comparison of simple potential functions for simulating liquid water.
        J. Chem. Phys. 1983; 79: 926
        • Beglov D.
        • Roux B.
        Finite representation of an infinite bulk system: solvent boundary potential for computer simulations.
        J. Chem. Phys. 1994; 100: 9050
        • Darden T.
        • York D.
        • Pedersen L.
        Particle mesh Ewald: an N log(N) method for Ewald sums in large systems.
        J. Chem. Phys. 1993; 98: 10089-10092
        • Essmann U.
        • Perera L.
        • Berkowitz M.L.
        • Darden T.
        • Lee H.
        • Pedersen L.G.
        A smooth particle mesh Ewald method.
        J. Chem. Phys. 1995; 103: 8577
        • Feller S.E.
        • Zhang Y.
        • Pastor R.W.
        • Brooks B.R.
        Constant pressure molecular dynamics simulation: the Langevin piston method.
        J. Chem. Phys. 1995; 103: 4613-4621
        • Martyna G.J.
        • Tobias D.J.
        • Klein M.L.
        Constant pressure molecular dynamics algorithms.
        J. Chem. Phys. 1994; 101: 4177
        • Kubo R.
        The fluctuation-dissipation theorem.
        Rep. Prog. Phys. 1966; 29: 255-284
        • Miyamoto S.
        • Kollman P.a.
        Settle: an analytical version of the SHAKE and RATTLE algorithm for rigid water models.
        J. Comput. Chem. 1992; 13: 952-962
        • Sugita Y.
        • Kitao A.
        • Okamoto Y.
        Multidimensional replica-exchange method for free-energy calculations.
        J. Chem. Phys. 2000; 113: 6042-6051
        • Fukunishi H.
        • Watanabe O.
        • Takada S.
        On the Hamiltonian replica exchange method for efficient sampling of biomolecular systems: application to protein structure prediction.
        J. Chem. Phys. 2002; 116: 9058-9067
        • Park S.
        • Kim T.
        • Im W.
        Transmembrane helix assembly by window exchange umbrella sampling.
        Phys. Rev. Lett. 2012; 108108102
        • Moradi M.
        • Tajkhorshid E.
        Computational recipe for efficient description of large-scale conformational changes in biomolecular systems.
        J. Chem. Theor. Comput. 2014; 10: 2866-2880
        • Torrie G.
        • Valleau J.
        Nonphysical sampling distributions in Monte Carlo free-energy estimation: umbrella sampling.
        J. Comput. Phys. 1977; 23: 187-199
        • Kumar S.
        • Rosenberg J.M.
        • Bouzida D.
        • Swendsen R.H.
        • Kollman P.A.
        The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method.
        J. Comput. Chem. 1992; 13: 1011-1021
        • Ferguson A.L.
        BayesWHAM: a bayesian approach for free energy estimation, reweighting, and uncertainty quantification in the weighted histogram analysis method.
        J. Comput. Chem. 2017; 38: 1583-1605
        • Metropolis N.
        • Rosenbluth A.W.
        • Rosenbluth M.N.
        • Teller A.H.
        • Teller E.
        Equation of state calculations by fast computing machines.
        J. Chem. Phys. 1953; 21: 1087-1092
        • Sindhikara D.J.
        • Emerson D.J.
        • Roitberg A.E.
        Exchange often and properly in replica exchange molecular dynamics.
        J. Chem. Theor. Comput. 2010; 6: 2804-2808
        • Darve E.
        • Pohorille A.
        Calculating free energies using average force.
        J. Chem. Phys. 2001; 115: 9169-9183
        • Comer J.
        • Gumbart J.C.
        • Hénin J.
        • Lelièvre T.
        • Pohorille A.
        • Chipot C.
        The adaptive biasing force method: everything you always wanted to know but were afraid to ask.
        J. Phys. Chem. B. 2015; 119: 1129-1151
        • Minoukadeh K.
        • Chipot C.
        • Lelièvre T.
        Potential of mean force calculations: a multiple-walker adaptive biasing force approach.
        J. Chem. Theor. Comput. 2010; 6: 1008-1017
        • Humphrey W.
        • Dalke A.
        • Schulten K.
        VMD: visual molecular dynamics.
        J. Mol. Graph. 1996; 14: 33-38
        • van der Walt S.
        • Colbert S.C.
        • Varoquaux G.
        The NumPy array: a structure for efficient numerical computation.
        Comput. Sci. Eng. 2011; 13: 22-30
        • Harris C.R.
        • Millman K.J.
        • van der Walt S.J.
        • Gommers R.
        • Virtanen P.
        • Cournapeau D.
        • et al.
        Array programming with NumPy.
        Nature. 2020; 585: 357-362
        • Stone J.E.
        • Vandivort K.L.
        • Schulten K.
        GPU-accelerated molecular visualization on petascale supercomputing platforms.
        in: Proceedings of the 8th International Workshop on Ultrascale Visualization. 2013https://doi.org/10.1145/2535571.2535595 (Denver, Colorado)
        • Hunter J.D.
        Matplotlib: a 2D graphics environment.
        Comput. Sci. Eng. 2007; 9: 90-95
        • Einstein A.
        Uber die von der molekularkinetischen Theorie der Wärmë geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen.
        Ann. Phys. 1905; 322: 549-560
        • Diamond J.M.
        • Katz Y.
        Interpretation of nonelectrolyte partition coefficients between dimyristoyl lecithin and water.
        J. Membr. Biol. 1974; 17: 121-154
        • Habeck M.
        Bayesian estimation of free energies from equilibrium simulations.
        Phys. Rev. Lett. 2012; 109100601
        • Sun L.
        • Bertelshofer F.
        • Greiner G.
        • Böckmann R.A.
        Characteristics of sucrose transport through the sucrose-specific porin ScrY studied by molecular dynamics simulations.
        Front. Bioeng. Biotechnol. 2016; 4https://doi.org/10.3389/fbioe.2016.00009
        • Golla V.K.
        • Prajapati J.D.
        • Joshi M.
        • Kleinekathöfer U.
        Exploration of free energy surfaces across a membrane channel using metadynamics and umbrella sampling.
        J. Chem. Theor. Comput. 2020; 16: 2751-2765
        • Hummer G.
        Position-dependent diffusion coefficients and free energies from Bayesian analysis of equilibrium and replica molecular dynamics simulations.
        New J. Phys. 2005; 7: 34
        • Gaalswyk K.
        • Awoonor-Williams E.
        • Rowley C.N.
        Generalized Langevin methods for calculating transmembrane diffusivity.
        J. Chem. Theor. Comput. 2016; 12: 5609-5619
        • Fujimoto K.
        • Nagai T.
        • Yamaguchi T.
        Momentum removal to obtain the position-dependent diffusion constant in constrained molecular dynamics simulation.
        J. Comput. Chem. 2021; 42: 2136-2144
        • Comer J.
        • Chipot C.
        • González-Nilo F.D.
        Calculating position-dependent diffusivity in biased molecular dynamics simulations.
        J. Chem. Theor. Comput. 2013; 9: 876-882