Formation of gaseous protein ions from aqueous ammonium acetate droplets during native electrospray: insights from mobile-proton molecular dynamics simulations

Abstract

The transition of proteins from solution into the gas phase during native electrospray ionization (ESI) continues to be an enigmatic topic. Molecular dynamics (MD) simulations are an important avenue for gaining mechanistic insights, but most previous MD studies did not consider the role of proton transfer (PT) events. Here, we overcome this limitation by implementing a droplet mobile-proton MD (dMPMD) protocol for aqueous ammonium acetate solutions. The droplet radius (up to 8.5 nm) and protein size (up to 55 kDa) are the largest studied to date, bringing the simulations closer to experimental conditions than any earlier modeling studies. Consistent with experiments, our data show that PT among NH₄⁺, acetate (Ac⁻), and protein sites produces NH₃ and acetic acid (HAc). NH₃ vaporizes rapidly, while HAc accumulation causes Ac⁻/HAc buffering at pH 4.76 ± 1. Droplet evaporation near the Rayleigh limit produces gaseous proteins with a net charge close to that of protein-sized water droplets, in line with the charged residue model (CRM). The CRM/Rayleigh framework applies regardless of amino acid sequence and residue basicity. For proteins with very few basic sites, charging involves carbonyl-trapped NH₄⁺ adducts. Switching to gas-phase mobile-proton MD (gMPMD) simulations, we track proteins after their release from droplets, culminating in native-like conformers with numerous salt bridges. Our data provide unprecedented insights into the ESI process, and the maturation of nascent protein ions to metastable structures that are detectable in experiments.

---

1 Introduction

Electrospray ionization (ESI) has become indispensable for mass spectrometry (MS) studies on proteins and other biomacromolecules.¹ “Native” ESI experiments aim to retain solution-like structures and interactions in the vacuum of the mass spectrometer.^2–10 Results from ion mobility spectrometry (IMS) and other techniques indicate that native-like protein conformations can indeed survive on experimental IMS/MS time scales,^11–16 provided that suitable conditions are used.^17,18 Some structural changes will nonetheless occur during the transition from solution into the gas phase. The exact nature of these changes remains a matter of debate.^19,20

The ESI process starts with solution in a high voltage capillary. Droplets emitted from a Taylor cone at the capillary outlet undergo repeated cycles of fission and evaporation. Throughout these events, the net number of charges on each droplet stays close to the Rayleigh limit^21–23

z_R = 8π/e × (ε₀γr³)^1/2 (1)

with r = droplet radius, γ = surface tension, e = 1.602 × 10⁻¹⁹ C, and ε₀ = 8.854 × 10⁻¹² F m⁻¹. Different mechanisms may be operational for the formation of desolvated protein ions from the final nanodroplets in the ESI plume. Previous work suggests that native ESI favors the charged residue mechanism (CRM), where proteins are released upon nanodroplet evaporation to dryness.^21,24,25 Alternatively, gaseous proteins may form via the ion evaporation mechanism (IEM) which involves the electrostatic ejection of ions from the nanodroplet surface.^26,27 The IEM dominates for small species such as metal cations, but even proteins can show IEM behavior under certain conditions.^25,28,29 Various other ion formation scenarios have been proposed as well.^30–35 Regardless of the mechanism, exposure to non-physiological pH or other unfavorable factors at any point during ESI can trigger protein unfolding and the dissociation of ligand interactions.^36–42

Bulk solution experiments promote the stability of native proteins by using neutral aqueous buffers with background electrolytes (e.g., 50 mM phosphate with 100 mM NaCl at pH 7). Unfortunately, such non-volatile solutes tend to associate with proteins during ESI, culminating in heterogenous adduct clusters instead of “clean” [M + zH]^z+ ions.^21,43 Several ESI-compatible (volatile) buffers and additives have been tested,^37,44–48 but most native ESI studies continue to rely on aqueous ammonium acetate (NH₄Ac).⁷ This salt decomposes via proton transfer (PT). The resulting products, ammonia and acetic acid, have a high evaporation propensity (NH₄⁺ + Ac⁻ → NH₃ + HAc).^21,43 Protein-bound NH₄⁺ can evaporate as NH₃, leaving behind a proton. Similarly, protein-bound Ac⁻ can evaporate as HAc after acquiring a proton.^4,21,49 Thus, NH₄Ac solutions favor the formation of [M + zH]^z+ ions without undesired adducts.²¹

Despite its widespread use, NH₄Ac is not an ideal additive for ESI. NH₄Ac produces a physiologically relevant near-neutral environment in water,^50,51 but NH₄Ac solutions do not have buffering capacity at pH 7.⁴³ Instead, they buffer ±1 pH units around the pK_a values of HAc (4.76) and NH₄⁺ (9.25). Thus, addition of acid or base to neutral NH₄Ac solution can cause large pH changes that destabilize native proteins.^38,40,43,46 Two factors contribute to acidification in positive ESI mode: (i) water oxidation in the emitter releases protons.⁵² (ii) For NH₄Ac-containing droplets, PT generates NH₃, a basic gas that evaporates rapidly. The other product (HAc) is a liquid acid that evaporates more slowly. Thus, HAc accumulation lowers the droplet pH until acetate buffering starts to take place, from the initial value of pH 7 down to 4.76 ± 1.⁴³ How exactly this mildly acidic droplet milieu affects proteins during the ESI process remains to be explored.^31,53

Molecular dynamics (MD) simulations of nanodroplets have become an important tool for probing ESI mechanisms.^35,54–65 Instead of focusing on [M + zH]^z+ ions that are most relevant in experiments, many earlier MD studies modeled analytes that were charged by metal adduction, e.g., [M + zNa]^z+. This approach sidesteps a limitation of classical MD force fields, i.e., their inability to capture the formation and dissociation of covalent bonds (such as PT events involved in [M + zH]^z+ formation). Over the past few years, mobile-proton MD (MPMD) methods have opened the door to simulations on ESI droplets^43,66,67 and gaseous protein ions^13,68,69 with PT.

We recently proposed an MPMD approach for modeling the behavior of NH₄Ac in ESI droplets (referred to as dMPMD, where the “d” stands for “droplet”).⁴³ For the first time, that work captured the PT-mediated formation of NH₃ and HAc. However, our study⁴³ only examined analyte-free ESI droplets. Our subsequent goal was to “extend these simulations to ESI droplets containing proteins […], to obtain a truly molecular understanding of native ESI”.⁴³ The present study accomplishes this goal.

While the current work was in progress, we became aware of a paper by Cordes and Gallagher.⁶⁷ Building on our recently developed methodology,^43,66 ref. 67 also modeled protein ion formation from NH₄Ac-containing aqueous droplets.⁶⁷ By pursuing a strategy complementary to that of ref. 67, the current work provides novel mechanistic insights into native ESI via a two-step strategy. Initially, we performed dMPMD simulations for modeling the release of protein ions from ESI droplets into the gas phase. Subsequently, the fate of these nascent ions was explored using gas-phase MPMD simulations (gMPMD, where the “g” stands for “gas phase”), to explore their structural evolution in vacuo. Our approach allows computational studies on droplets and proteins in an unprecedented size range. For the first time, these simulations capture the complete sequence of ESI events, from proteins in droplets to metastable gas-phase structures that can be compared to IMS experiments, with inclusion of PT events along the entire pathway.

2 Methods

ESI-IMS/MS experimental settings and MD details can be found in the SI.

2.1 Droplet simulations: general aspects

Like previous ESI simulations,^58,65 we used the CHARMM36 force field.⁷⁰ The application of pseudo-periodic boundary conditions (for droplets and for gas phase proteins, see Section 2.4) ensured an environment equivalent to vacuum boundary conditions,⁵⁵ thereby eliminating artifacts caused by finite cutoffs for nonbonded interactions.⁷¹ Protein titratable sites were initially set for pH 5.4,⁴³ i.e., nt⁺ (N-terminus), R⁺, K⁺, H⁺, D⁻, E⁻, and ct⁻ (C-terminus) (Fig. S1). Crystal structures of ubiquitin (1UBQ), lysozyme (1AKI), and pepsin (4PEP) served as starting points. Missing residues in TTR file 3GRG necessitated the use of AlphaFold 3 to generate starting structures (Fig. S3).⁷²

Droplets contained water, NH₄⁺, Ac⁻, and HAc, with a protein initially centered. Parameters for non-protein solutes were taken from ref. 43. Molalities were used to ensure consistent NH₄Ac concentrations in the non-protein part of the droplets. Ubiquitin droplets with an initial radius of r₀ = 3 nm had a NH₄Ac molality of 2.3 m, corresponding to a molarity of 1.6 M when considering the entire droplet volume. For 3 nm lysozyme droplets, we used 2.3 m (1.5 M). Molalities of 6 nm droplets were 1.15 m (1.00 M, 0.99 M, 0.96 M, and 0.93 M for ubiquitin, lysozyme, pepsin, and TTR). Ubiquitin runs with r₀ = 8.5 nm employed 1.05 m (0.92 M). These NH₄Ac concentrations are higher than those of bulk ESI solution, where [NH₄Ac] is typically 10–100 mM.⁷ However, the droplets simulated here represent a regime where water evaporation has increased solute concentrations.²¹ The [NH₄Ac] of 1.6 M for r₀ = 3 nm matches the conditions of our earlier work.⁴³ NH₄Ac concentrations for larger droplets were chosen such that [NH₄Ac] approached 1.6 M once the droplets had shrunk to 3 nm. Additional HAc were included to ensure an initial [Ac⁻]/[HAc] ratio consistent with the Henderson-Hasselbalch equation (eqn (2)) for pH 5.4.⁴³

Additional NH₄⁺ were inserted to bring the droplet charge to the Rayleigh limit (eqn (1)) with γ = 0.05891 N m⁻¹. The resulting z_R values were 18+ for r₀ = 3 nm, 52+ (6 nm), and 88+ (8.5 nm). Negatively charged droplets were modeled with excess Ac⁻ to implement the droplet charge, while adjusting the initial number of HAc to ensure a pH of 5.4 (eqn (2)).

Droplets underwent energy minimization, followed by 500 ps of equilibration during which the temperature was raised to 370 K. After maintaining 370 K for much of the run, the temperature was raised to 450 K for promoting the final evaporation steps. The 370 K/450 K periods were 50 ns/25 ns for r₀ = 3 nm, 112.5 ns/37.5 ns (6 nm), and 137.5 ns/37.5 ns (8.5 nm). Although typical ESI interfaces employ droplet heating, the MD temperatures used here are likely higher than those in experiments, where droplet temperatures have been estimated to be ∼320 K.⁷³ Unfortunately, droplet runs at T ≤ 320 K yield slow evaporation rates that preclude studies of ESI events.⁵⁵ This is especially true for TIP4P/2005 water which evaporates more slowly than other models, keeping in mind that those other models have unrealistically low surface tension.⁷⁴ Using TIP4P/2005 water at 370 K/450 K overcomes these problems, as shown previously.^60,61,67,71 A possible concern with the temperature regime in our simulations is the occurrence of thermal protein unfolding. However, the data discussed below (as well as earlier experiments⁷⁵ and simulations^60,61,67,71) demonstrate that the short time scale of ESI precludes thermal unfolding under the conditions used here. Droplet volume was determined as V_droplet = [mass of water + solutes + protein (kg)]/[1000 kg m⁻³], and droplet radii for eqn (1) were estimated using V_droplet = 4/3πr³.

2.2 Droplet pH

The pH of evaporating droplets was determined from their Ac⁻/HAc ratio (eqn (2)).⁴³ The small droplet size studied here implies that the pH-derived number of free H₃O⁺ (N_H3O⁺) tends to be less than 1. Analogous situations are encountered in bioenergetics, when studying sub-cellular compartments.^76,77 The definition pH = –log[H₃O⁺] remains valid in such small-volume systems, provided that N_H3O⁺ is interpreted as a temporal average.^76,77 Specifically, one can consider the average number 〈N_H3O⁺〉 = V_droplet × 10^−pH × N_A for droplets undergoing rapid Poissonian fluctuations between instances where they contain 0, 1, … H₃O⁺.^76,77 For example, a droplet containing zero H₃O⁺ (90% of the time) and 1H₃O⁺ (10% of the time) would have 〈N_H3O⁺〉 = 0.1. The situation is even more striking for the pH range considered here; e.g., a droplet with radius 3 nm at pH 3.76 (the lower end of Ac⁻/HAc buffering) has 〈N_H3O⁺〉 ≈ 0.01. The fleeting nature of H₃O⁺ under these conditions prompted us to model ESI droplets without explicitly including H₃O⁺.

2.3 Droplet MPMD (dMPMD) simulations

Simulations used trajectory stitching, where runs were broken down into 250 ps segments. Evaporated species (>4 × r₀ from the center) were removed after each segment. The next segment then started with the remaining coordinates and velocities. Repeat runs used different initial coordinates and velocities. After each segment, the algorithm allowed for PT from protonated donors to Ac⁻, and from NH₄⁺ to deprotonated acceptors (Fig. 1 and Fig. S4). Possible donors were NH₄⁺ (any HN atom), nt⁺ (any HN atom), R⁺ (HNδ), K⁺ (any HN atom), H⁺ (HNδ1 or HNε2), D⁰, E⁰, and ct⁰. Possible acceptors were Ac⁻, nt⁰, R⁰ (Nδ), K⁰ (Nε), H⁰ (deprotonated Nδ1 or Nε2). Also, any oxygen in D⁻, E⁻, and ct⁻ was a possible acceptor. PT among side chains (or side chains/termini) was not considered because it generally exceeded the maximum allowed distance.

Fig. 1 Proton donors and acceptors considered for the dMPMD algorithm of this work. Arrows indicate possible PT reactions. Snapshots of all PT events are exemplified in Fig. S4.

Donor/acceptor pairs were identified as PT candidates if they were H-bonded to one another, with a distance less than a specific cutoff (PT_dist_max, see below). Unless noted otherwise, PT was executed with a probability PT_prob = 0.1 for each candidate pair. NH₄⁺/NH₃ and Ac⁻/HAc conversion was performed as described⁴³ (Fig. S4A, top panel). Protonation/deprotonation of side chains and termini was performed using GROMACS pdb2gmx. Rapid NH₃ evaporation from the droplet justifies modeling PT from NH₄⁺ to Ac⁻ as an irreversible step.⁴³ Protein sites were able to switch back and forth between protonation states, e.g. K⁺ deprotonation by Ac⁻, with subsequent K⁰ protonation by NH₄⁺. For determining the maximum distance for which PT was permitted (PT_dist_max), we generated H-bond distance histograms for all possible donor/acceptor pairs (Fig. 1) using the strategy of Fig. S5. Histogram maxima were used as PT_dist_max for each pair (Fig. S6).

2.4 Gas-phase MPMD (gMPMD) simulations

Desolvated protein ions released from ESI droplets were subjected to gMPMD.^13,78 The temperature was set to 300 K to mimic non-activating IMS/MS conditions.^{2–7,11–18} To ensure consistency with previous gMPMD work,^13,78 the OPLS-AA force field⁷⁹ was used. 200 ns trajectories were stitched together from 1 ns segments with fixed proton configurations. After each segment, protonated sites (nt⁺, R⁺, K⁺, H⁺, E⁰, D⁰, ct⁰) were allowed to undergo PT to deprotonated sites (nt⁰, R⁰, K⁰, H⁰, E⁻, D⁻, ct⁻). The next segment then continued from the updated proton configuration. Proton redistribution was governed by minimization of the energy term

E_gMPMD = E_Coul + E_PA (3)

E_Coul captures intramolecular charge solvation via electrostatic interactions among all atomswhere q_i and q_j are atomic charges, and r_ij are distances. Protonation of a titratable site lowers E_MPMD by an amount equal to its intrinsic proton affinity (PA).⁸⁰ For N sites, this corresponds towith δ_k = 1, 0 for protonated and deprotonated sites, respectively. Gly was added to the lysozyme nt and the TTR ct to avoid the presence of two titratable sites in one residue.⁷⁸ gMPMD runs started from droplet-generated protein structures and charge states (z). To account for declustering during ion sampling,²¹ residual adducts were removed without altering z (e.g., NH₄⁺ was replaced by a proton). From a random protonation configuration with the appropriate z, proton redistribution was then performed by minimizing E_gMPMD (eqn (3)), before commencing the gMPMD runs.

2.5 Comparison of algorithms

Compared to the dMPMD/gMPMD approach of ref. 67, the strategy used here is conceptually simpler and has lower computational cost, allowing simulations in an unprecedented size range (see Section 3.4 and 3.5). Key differences compared to ref. 67 are as follows:

(i) Ac⁻/HAc buffering stabilizes the droplet pH at 4.76 ± 1 in native ESI experiments and in our simulations (eqn (2)).⁴³ The sparsity of H₃O⁺ under these conditions prompted us to model droplets without explicitly including H₃O⁺ (see Section 2.2). This strategy eliminated the need to simulate H₃O⁺ Grotthuss diffusion,^66,67 thereby streamlining the code and facilitating its application to large systems. In contrast, ref. 67 modeled ESI droplets with a high H₃O⁺ abundance that is difficult to reconcile with a Ac⁻/HAc-buffered milieu.⁴³ It cannot be ruled out that in ESI experiments a small fraction of protein protonation events involve transiently formed H₃O⁺. However, for any H₃O⁺ that might be present, Ac⁻ is a much more likely proton acceptor (e.g., ubiquitin only has ≤12 vacant protonation sites at pH 5.4 [Fig. S1], while a r₀ = 6 nm droplet contains 544 Ac⁻). Thus, the direct participation of H₃O⁺ in protein protonation will be minimal, providing additional justification for not explicitly considering H₃O⁺ in our MD strategy.

(ii) Ref. 67 executed PT on the basis of ΔG° estimates. Their eqn (2) predicts that NH₄⁺ + Ac⁻ → NH₃ + HAc in droplets has ΔG° = +26 kJ mol⁻¹. This unfavorable value is at odds with the fact that NH₃ and HAc readily form in evaporating NH₄Ac solutions,⁴³ illustrating that bulk thermodynamics may be unsuitable for predicting the behavior of non-equilibrium systems with liquid/vapor interfaces. Similarly, ΔG° values estimated from gas-phase basicities (GPBs) in ref. 67 did not consider electrostatic charge solvation effects which are key for determining PT favorability.^13,53,78,80 The use of the empirical probability parameter PT_prob in our work sidesteps these issues (see Section 2.3).

3 Results and discussion

3.1 Assessing the viability of dMPMD simulations

The goal of this work was to design an efficient strategy for modeling native ESI, i.e., the formation of gaseous protein ions from aqueous NH₄Ac droplets, with inclusion of PT. This was achieved by extending our earlier dMPMD approach, where proteins were not considered.⁴³ In addition to NH₄⁺ and Ac⁻, all protein titratable sites were allowed to participate in PT (Fig. 1). As a first step, we verified the fidelity of our computational strategy by applying it to relatively small droplets (r₀ = 3 nm), representing a regime close to the end of the fission/evaporation cascade in the ESI plume.²¹ Also, this size range matches typical dimensions used in earlier MD work.^35,54–63

Ubiquitin (8656 Da, pI 6.6) is a widely studied model protein.^11,67,78,81 Fig. 2A illustrates trajectory snapshots in positive ion mode. Water evaporation caused droplet shrinkage. PT-generated NH₃ and HAc evaporated as well. Charge loss from the evaporating droplets took place via NH₄⁺ IEM ejection.^26,27 These events culminated in 6+ gaseous ubiquitin. The observed droplet evaporation to dryness implies that this protein ion was formed via the CRM.^21,24

Fig. 2 Ubiquitin native ESI (dMPMD) simulations in r₀ = 3 nm aqueous NH₄Ac droplets, using positive ESI mode. (A) Trajectory snapshots. NH₄⁺, Ac⁻, NH₃, and HAc are shown in spacefill representation, water is shown as lines that are barely discernible, the protein is colored green. All other panels show data from 20 repeat runs. (B) Number of H₂O in droplet, (C) net droplet charge, (D) number of solutes in droplet, (E) cumulative number of PT events from NH₄⁺ to Ac⁻, and PT events involving the protein, (F) droplet pH, (G) relative droplet charge. Panels F and G only show the time range down to ∼50 H₂O, because the corresponding eqn (1) and (2) refer to aqueous systems. Also included in panel G is the average z/z_R in this 20 ns window. (H) Protein adducts at the end of droplet evaporation. (I) Simulated protein ESI charge state distribution (including adducts).

Fig. 2 includes data from 20 repeat runs, all of which followed a pattern similar to Fig. 2A. Most water evaporated within 25 ns (Fig. 2B), concomitant with a stepwise decrease in droplet charge via NH₄⁺ IEM events (Fig. 2C). PT produced NH₃ and HAc (Fig. 2D and E). The high volatility of NH₃ translated in a low NH₃ steady-state concentration (Fig. 2D). In contrast, the lower volatility of HAc caused accumulation of this acid, resulting in a pH decrease from 5.4 to ca. 4.5 (Fig. 2D, F and eqn (2)). The droplet charge remained around 80% of the Rayleigh limit (Fig. 2G and eqn (1)). Some proteins retained one or two solutes at the end of the 75 ns simulations, with binding of a single NH₄⁺ as the most common outcome (Fig. 2H). The ESI charge state distribution, obtained by tallying all charges at the end of the 20 runs (side chains, termini, and adducts) ranged from 5+ to 7+, consistent with experiments^67,78 (Fig. 2I).

PT of donor-acceptor candidate pairs was executed with a probability PT_prob (see Section 2.3). This probability represents an empirical parameter rather than a physically derived property. The value of PT_prob = 0.1 used for Fig. 2 matches the settings employed earlier for protein-free droplets.⁴³ To test the effects of PT_prob, we performed ubiquitin simulations under different PT_prob regimes (Fig. S7). The different settings affected how the droplet composition changed over time. However, all simulations ultimately produced similar protein charge states (Fig. S7, bottom), covering the range of 5+ to 7+ that matches experimental spectra.^67,78 In other words, the simulation outcomes were largely insensitive to the choice of PT_prob, prompting us to retain PT_prob = 0.1 for all data discussed below.

Simulations were also conducted on lysozyme (14 305 Da, pI 10.5), another common model protein.^67,78 ESI events for lysozyme resembled those for ubiquitin (Fig. S8). dMPMD-generated lysozyme charge states ranged from 7+ to 9+, matching experimental data.^67,78

3.2 Positive and negative ion mode ESI

While most ESI studies employ positive ion mode, experiments can also be performed in negative ion mode where protein anions are released from negatively charged droplets.^67,78,82 After verifying the viability of our algorithm in positive ion mode (Fig. 2 and Fig. S8), we simulated r₀ = 3 nm droplets in negative ion mode. The net droplet charge was implemented by excess Ac⁻. The ensuing events were analogous to those in positive ion mode, i.e., evaporating droplets underwent Ac⁻ IEM events, yielding CRM charge states that matched the experimentally observed range (5- to 6- for ubiquitin, 6- to 8- for lysozyme, Fig. S9 and S10).^67,78

Native ESI experiments in negative ion mode tend to produce slightly lower charge states |z| than in positive ion mode.^67,78,82,83 This effect was reproduced in our simulations, where average z values were 5.5- vs. 6.1+ for ubiquitin, and 7.0- vs. 8.1+ for lysozyme (Fig. 2 and Fig. S8–S10). Our simulations reveal that this disparity is rooted in the ejection behavior of excess charge carriers, NH₄⁺ vs. Ac⁻, illustrated for ubiquitin in Fig. S11. Charge ejection always occurred with a solvation shell. For droplets that still contained abundant water, ejected NH₄⁺ and Ac⁻ were mainly solvated by H₂O (Fig. S11C and J). As the droplets gradually dried out, solvation by HAc became more prevalent (Fig. S11G and N). The absolute system charge in positive and negative ion mode was indistinguishable at t ≈ 30 ns (Fig. S11B), implying that “wet” droplets do not hold the key for explaining the charge disparity between positive and negative ion mode. Instead, differences manifested themselves when the systems were almost anhydrous. Charge loss in positive ion mode became rare at t > 30 ns, with only a single ejection in 20 simulations (Fig. S11C–G). Negative ion mode exhibited a higher charge loss propensity at this late stage, with 14 ejections for t > 30 ns (Fig. S11J–N). Thus, our data show that toward the end of droplet evaporation, Ac⁻-mediated charge loss in negative ion mode is more facile than NH₄⁺-mediated charge loss in positive ion mode. This trend was also seen for lysozyme (Fig. S12).

In summary, the dMPMD strategy developed here is suitable for modeling native protein ESI in positive and negative ion mode. All r₀ = 3 nm droplet runs showed protein ion production via the CRM, accompanied by IEM charge ejection (solvated NH₄⁺ or Ac⁻ in positive or negative ion mode, respectively). Simulated protein charge states were in agreement with experiments.⁷⁸ Except for Fig. S9–S12, all data in this work represent positive ion mode.

3.3 dMPMD simulations on r₀ = 6 nm droplets

The droplets discussed above with r₀ = 3 nm represent a late stage of the ESI process.²¹ Many earlier studies were restricted to this size range, as dictated by computational cost.^{35,54–63,65} The streamlined approach used here facilitates investigations on much larger droplets, offering the opportunity to capture a more extensive segment of the ESI process. Larger droplets are also required for modeling high MW proteins that exceed the size of r₀ = 3 nm droplets.

dMPMD simulations on r₀ = 6 nm droplets were performed for ubiquitin and lysozyme, as well as pepsin which represents a large monomeric protein (35 kDa, pI ≈ 3, protein radius ≈ 3 nm). We also modeled droplets containing transthyretin (TTR), a 55 kDa noncovalent tetramer (pI ≈ 5, protein radius ≈3.5 nm). Results for TTR are illustrated in Fig. 3. Data for the other three proteins are summarized in Fig. S13–S15. All r₀ = 6 nm trajectories showed droplet evaporation to dryness, following a sequence of events similar to that seen for the r₀ = 3 nm droplets in Fig. 2 and Fig. S8–S10. Simulated charge state distributions for ubiquitin were virtually identical for 3 nm and 6 nm droplets (Fig. 2I and Fig. S13I). The same was true for lysozyme (Fig. S8I and Fig. S14I). Overall, our simulations reveal that the CRM is operative for proteins across a wide size range (we tested 8.5 kDa to 55 kDa), with different acid/base properties (we examined pI 3 to 10.5), and for droplet radii up to at least 6 nm. TTR represents the largest protein to date that has been subjected to ESI simulations (Fig. 3).

Fig. 3 TTR native ESI (dMPMD) simulations in r₀ = 6 nm aqueous NH₄Ac droplets, using positive ESI mode. (A) Representative trajectory snapshots. Panels B-I summarize data from 20 replicates. For additional details, see the caption of Fig. 2.

Simulated ESI charge state distributions for ubiquitin, lysozyme, and pepsin from r₀ = 6 nm droplets were almost identical to the results of native ESI-MS experiments (Fig. 4A–F). Published TTR spectra show some variability, reflecting the TTR propensity to undergo chemical modifications (Fig. S2), and the capability of the flexible N-termini to adopt different conformations (Fig. S3). Experimental TTR charge states in the literature cover the range of 13+ to 17+.^17,84,85 Our experiments yielded 15+ to 17+ (Fig. 4G). The simulated TTR charge states of 14+ to 16+ are consistent with the experimentally observed range (Fig. 4H). The overall agreement between experimental spectra and simulation results in Fig. 4 is remarkable, attesting to the robustness of our dMPMD approach.

Fig. 4 Native ESI mass spectra compared to dMPMD-generated protein charge states for r₀ = 6 nm droplets. (A) and (B) Ubiquitin, (C) and (D) lysozyme, (E) and (F) pepsin, and (G) and (H) TTR. Each simulated spectrum is based on 20 replicates. Dashed red lines and red numbers indicate ESI charge states predicted as 0.85 × z_R, where z_R is the Rayleigh charge of a protein-sized water droplet (eqn (1)), using radii determined from the protein mass assuming spherical shape and a 1 g cm⁻³ density.²⁴ Shown along the right are representative protein structures after 150 ns of droplet evaporation.

3.4 Residual adducts

Native ESI experiments on small proteins often produce ions that are adduct-free (Fig. 4A and C). This behavior arises from the PT events in Fig. 1, that convert non-volatile NH₄⁺ and Ac⁻ into volatile NH₃ and HAc.^21,43 dMPMD-generated protein ions carried slightly higher adduct levels than in experiments, with NH₄⁺ as the most prevalent charged adduct. Ubiquitin and lysozyme mainly formed [M + (z-1)H + NH₄]^z+. NH₄⁺ binding resulted from “carbonyl trapping” by backbone and side chain C=O groups (Fig. 5A).

Fig. 5 Carbonyl trapping of NH₄⁺. (A) NH₄⁺ binding to ubiquitin under the conditions of Fig. 2. Red labels indicate backbone (bold) and side chain (italicized) oxygens. (B) Pepsin 11+ produced in the dMPMD runs of Fig. S15, with 14 bound NH₄⁺ ions. (C) and (D) Close-up of NH₄⁺ binding to pepsin. (E) and (F) MS/MS data acquired with gentle activation (trap CV 2 V, panel E) and with moderate activation (15 V, panel F). Theoretical isotope distributions of [pepsin + (11-n)H + nNH₄]¹¹⁺ are shown for n = 0 (red), and n = 1, …, 5 (blue). Panel F also shows water loss, a common event during collisional activation of S/T-containing proteins (the S/T content of pepsin is 21%).⁸⁷

Carbonyl-trapped NH₄⁺ was particularly abundant in pepsin simulations, with occasional participation of other binding sites such as Ser-OH (Fig. 5B–D and Fig. S15H). Consistent with their high abundance in our MD data, NH₄⁺ adducts were observed in native ESI mass spectra of pepsin (Fig. 5E and Fig. S16). Collisional activation caused a gradual shift to lower m/z (Fig. 5E and F). The process did not change the protein charge, implying NH₃ loss rather than NH₄⁺ ejection (Fig. S17). This interpretation is consistent with earlier experiments.^4,49 The MS/MS data of Fig. 5 suggest that collisional heating triggers PT from residual NH₄⁺ to acceptor sites on the protein (acceptor + NH₄⁺ → acceptor-H⁺ + NH₃).²¹ The scarcity of canonical basic acceptors in pepsin (a total of 5 sites, Fig. S1) implies the involvement of relatively unfavorable acceptors after collisional activation, such as side chains of Asn, Gln, Pro, and possibly even backbone carbonyls.⁸⁶ These unfavorable proton acceptors likely also participate under denaturing conditions that produce [pepsin + zH]^z+ ions without NH₄⁺ adducts (Fig. S2D).

Native ESI mass spectra of large protein complexes tend to show peak broadening, an effect often attributed to “salt adduction”.⁷ Adduct-mediated broadening is exemplified by the TTR experimental data in Fig. 4G. Our dMPMD results show that residual HAc is a main contributor to peak broadening (Fig. 3H). Some HAc adducts were buried in internal cavities of TTR, illustrating why native ESI can produce adducted ions even in the absence of salt contaminants (Fig. 3A, t = 150 ns). Collisional activation during ion sampling promotes partial loss of these adduct (Fig. S18).

3.5 Pushing the MD size range to r₀ = 8.5 nm droplets

Ideally, ESI simulations would start with droplets in the r₀ ≈ 100 nm range, as this size is produced at the Taylor cone in typical experiments.²¹ In addition to capturing analyte ion formation, such simulations would be expected to show droplet fission events that precede formation of the final nanodroplets. Sadly, r₀ ≈ 100 nm droplets are far out of range for existing MD methods. The minimum radius required for droplet fission has been estimated to be ca. 8.5 nm, whereas charge loss for smaller droplets is limited to IEM events²⁶ (as seen here for r₀ = 3 nm and 6 nm).

To access a regime that covers a wider range of behaviors, we performed dMPMD on ubiquitin-containing droplets with r₀ = 8.5 nm (∼76 000 H₂O). To our knowledge, these are the largest ESI droplets simulated using MD techniques (Fig. 6A). Three runs were performed, each of which required ∼1 month of wall clock time. Each run produced a gaseous protein via the CRM (Fig. S19), similar to the smaller droplets discussed earlier (Fig. 2 and Fig. S13). Protein charge states were 5+ (twice) and 7+, consistent with the experimentally observed range (Fig. 4).

Fig. 6 (A) Initial droplet size and water content for the different dMPMD simulations of this work. Ubiquitin is shown for size comparison. All other panels illustrate evaporation and jet fission of a r₀ = 8.5 nm droplet containing NH₄Ac and ubiquitin in positive ion mode. (B) and (C) Absolute and relative droplet charge; arrows indicate the fission time point. (D) and (E) dMPMD snapshots marking the onset and the end of jet fission. Panel E shows the composition and charge of five progeny droplets (highlighted by circles). Full ESI trajectories are shown in Fig. S19.

The shrinking r₀ = 8.5 nm droplets retained a spherical shape throughout most of their life cycle (Fig. S19A). NH₄⁺ and Ac⁻ were distributed uniformly in the aqueous phase (Fig. S20A–C). The NH₃ radial distribution showed a local maximum at the droplet surface, representing PT products that had diffused to the droplet periphery and that were about to evaporate (Fig. S20D). Highly pronounced surface accumulation was seen for HAc (Fig. S20E). The high surface affinity reflects the less favorable water solvation of HAc compared to Ac⁻. HAc enrichment at the liquid/vapor interface is consistent with the ESI behavior of other partially nonpolar species.^88,89 Most of the HAc molecules at the surface had their carboxylic acid group in contact with water, while the desolvated methyl group protruded into the vapor phase (Fig. S20H). The protein tended to reside relatively close to the surface as well, instead of retaining its initial centered position in the droplet (Fig. S20F).

Strikingly, each of the large droplets underwent a jet fission event, causing a sudden reduction in overall charge (Δz = 9 in Fig. 6B and C). Fission started with elongation of the initially round droplet (Fig. 6A, t = 0 in Fig. S19A), followed by formation of a Taylor cone-like protrusion that released a burst of progeny droplets from its tip (Fig. 6D and E). Several progenies were multiply charged; the largest one in Fig. 6E had a 3+ charge and contained ∼200 molecules. This progeny size is in contrast to the IEM ejection of singly charged ions (typically NH₄⁺ with ∼10 H₂O) that represents a competing charge loss mechanism (Fig. 2A and 3A and Fig. S19A). None of the progeny droplets contained a protein under the conditions studied here, although we cannot rule out the feasibility of such protein IEM or “intermediate regime ESI” events.^25,31

The MD data in Fig. 6 resemble experimentally observed jet fission events.^21–23 Only very few past MD studies ventured into a droplet size regime where such events become feasible.^90,91 The unprecedented size used here for the first time produced trajectories that captured the entire range of droplet behaviors, i.e., jet fission, evaporation to dryness with numerous IEM events, culminating in a CRM-generated gaseous protein ion (Fig. 6 and Fig. S19, S20).

Simulation data supporting the CRM were previously obtained for much smaller droplets.⁶⁵ A weakness of those earlier studies is the possibility that their outcomes might have been biased by the initial droplet size used. It is gratifying that trajectories consistent with the CRM were also obtained for the very large droplets simulated here, confirming that earlier mechanistic insights were not tainted by size-related artifacts.

3.6 The CRM/Rayleigh framework

Our dMPMD data reveal that under native ESI conditions each protein remains trapped inside a shrinking droplet, up to the point where the solvent has evaporated to dryness. This scenario matches the definition of the CRM.⁹² Evaporating ESI droplets are governed by two competing trends: (i) H-bonds and other cohesive interactions favor a near-spherical droplet shape that maximizes intermolecular contacts by minimizing the surface-to-volume ratio. (ii) Electrostatic repulsion destabilizes the droplet, promoting charge loss in the form of droplet fission and IEM events. The competing factors (i) and (ii) cause the evaporating droplets to stay in a regime where the net droplet charge is at ca. 85% of the Rayleigh limit (0.85 is the average z/z_R from all r₀ = 6 nm runs). Evaporation ends when only the dry protein remains, at which point the charge of the vanishing droplet becomes the charge of the nascent protein ion. Droplet shrinkage at ∼0.85 z_R, implies that the protein charge will be z ≈ 0.85 z_R of a protein-sized droplet. The predictions of this CRM/Rayleigh framework are in good agreement with simulated and experimental ESI charge states (Fig. 4, dashed lines), supporting earlier proposals.^21,24,49,93

Fig. 4 demonstrates that the CRM/Rayleigh framework applies to proteins that have vastly different sizes and basicities. An interesting scenario arises for proteins that possess very few basic sites. Specifically, pepsin contains only five canonical basic moieties (1nt, 2R, 1K, 1H, Fig. S1). From this amino acid composition, one might expect that pepsin should be unable to carry a charge greater than 5+. However, the CRM/Rayleigh framework predicts a native ESI charge of ∼11+, which is in agreement with experiments and simulations (Fig. 4E and F). Our data show that native ESI charging beyond 5+ in pepsin is mediated by carbonyl-trapped NH₄⁺ ions (Fig. 5 and Fig. S15H), as previously suggested on the basis of experimental data.⁴⁹

It has been proposed that residue basicity plays a key role for protein charging in native ESI.⁶⁷ This assertion contradicts experiments that demonstrated major disparities between analyte basicity and ESI charge states.⁹⁴ Instead, the CRM/Rayleigh framework implies that protein size is the single most important determinant of native ESI charge states.^24,49,93,94 The expected size dependence was reproduced in our simulations for proteins representing a wide range of basicities, demonstrating that basicity is not a key determinant of native ESI protein charge states.^24,49,93,94

3.7 From nascent ions to metastable gas-phase structures

After their release from ESI droplets, protein ions adapt to the solvent-free environment.^11,20 The net charge of each ion is fixed once it is in vacuo, but spatial protonation patterns continue to evolve because protons are highly mobile in the gas phase.^68,95–98 Many earlier MD studies have explored gaseous proteins,^{11,53,68,69,99–104} but most of them faced two limitations. (i) They started from X-ray or NMR structures, although droplet-generated structures would be a more realistic starting point. (ii) Mobile protons were not considered. Here we address both of these issues. (i) Gaseous proteins were modeled as a continuation of droplet simulations. (ii) We allowed for PT among all titratable sites, using a gas-phase mobile proton MD (gMPMD) algorithm where PT is governed by proton affinities and electrostatics (see Section 2.4).^13,78

Following their CRM production from r₀ = 6 nm droplets (Section 3.3), ubiquitin, lysozyme, and TTR in the most abundant charge states were modeled by gMPMD. Fig. 7 illustrates data for TTR. The radius of gyration (R_g, reflecting the overall protein size) remained close to the initial solution value throughout the droplet stage of the simulations (0–150 ns). During the subsequent 200 ns in the gas phase, TTR underwent a slight compaction, culminating in a R_g 6.5% lower than for the solution structure (Fig. 7A). Similar compaction events have been reported for other complexes.^100,101,103 Structural changes throughout the droplet and gas-phase simulations resulted in a ∼1 nm RMSD relative to the solution structure (Fig. 7B). Fig. 7C puts this RMSD into context, highlighting that most aspects of the solution architecture were retained in the gas phase, albeit with some distortion of β sheets and loops.

Fig. 7 Structural evolution of [TTR + 15H]¹⁵⁺ during CRM release from r₀ = 6 nm droplets (0–150 ns), and subsequently in the gas phase (150–350 ns). (A) Radius of gyration R_g. The initial solution R_g from Fig. S3 is shown as well. (B) Backbone RMSD relative to the initial solution structure. (C) Initial solution structure and representative final gas-phase structure. (D) Typical final gas-phase structure with highlighted positive (blue) and negative sites (red). Labels are shown only for one subunit. (E) Initial charge pattern in the droplet, and at the start/end (F) and (G) of gas-phase simulations. Error bars illustrate the variability of charge patterns across all subunits in the 10/20 runs highlighted here; the remaining runs generated charge states other than 15+.

The final gaseous TTR structures exhibited an extensive salt bridge network, involving protonated basic sites and R-COO⁻ moieties (Fig. 7D). The presence of these R-COO⁻ moieties in [M + zH]^z+ ions may seem counterintuitive, prompting many earlier MD studies to employ all-positive protonation patterns. However, it is now well established that intramolecular charge solvation and nearby positive charges can stabilize R-COO⁻ groups (eqn (4)).^53,105,106 The zwitterionic motifs formed in our gMPMD simulations are consistent with those data (Fig. 7D). Starting from the initial pH 5.4 droplet pattern (Fig. 7E), the zwitterionic character of TTR initially diminished until the droplet stage of the simulations had terminated (Fig. 7F). During the subsequent 200 ns in the gas phase, intramolecular PT generated additional positively charged sites and R-COO⁻ groups. All of the latter were involved in salt bridges (Fig. 7G and Fig. S21).

gMPMD simulations on ubiquitin and lysozyme yielded data similar to those for TTR, with numerous salt bridges, and with retention of native-like structures (Fig. S21–S23). A few ubiquitin and lysozyme runs generated slightly unfolded conformers, indicating that the gas phase resilience of these low MW proteins is somewhat lower than that of TTR, likely reflecting their larger surface-to-volume ratio. The propensity of ubiquitin to produce a small fraction of semi-unfolded conformers has been noted previously in IMS experiments,^11,78 although most ions retained compact structures (consistent with our simulations, Fig. S22).

In the case of pepsin, the scarcity of basic residues and the overabundance of acidic sites called for a slightly different simulation strategy. Consistent with experiments,⁴⁹ dMPMD-generated ions showed complete protonation of acidic sites. Protonation of the five basic sites, along with six NH₄⁺ adducts, resulted in a 11+ overall charge. Complete protonation of carboxylates implies the absence of salt bridges and mobile protons. Gas-phase simulations on pepsin were therefore performed with static protons. Pepsin retained a native-like but slightly collapsed gas-phase structure (5.9% lower R_g than the X-ray structure, Fig. S24).

The gMPMD structures of ubiquitin and lysozyme were in excellent agreement with IMS experiments (Fig. S25A–D). For pepsin and TTR, experimental collision cross sections were slightly lower than for the gMPMD structures, indicating that gas-phase collapse of these proteins did not fully go to completion on the 350 ns time scale of our simulations (Fig. S25E–H).

4 Conclusions

The mechanism(s) of protein ion formation during ESI, and the extent to which these ions retain solution-like structures have remained controversial for many years.^{19,20,25,31,53,67} Our simulations offer unprecedented insights into these events, focusing on NH₄Ac droplets that represent native ESI conditions, with inclusion of PT reactions. Our data capture the entire sequence of events, from proteins in droplets to nascent desolvated ions to IMS/MS-detectable gas-phase conformers. The droplets modeled here have radii up to 8.5 nm, the largest size simulated to date. Similarly, TTR (55 kDa) is the largest protein for which ESI simulations have been performed.

The dMPMD approach used here captures the experimental finding⁴³ that PT converts NH₄Ac solution into an Ac⁻/HAc buffer that stabilizes the droplet pH at 4.76 ± 1 (see, e.g., Fig. 2F and 3F). Our data demonstrate that this mildly acidic milieu allows the retention of folded protein conformations on the short (ns to μs)²¹ time scale of nanodroplet evaporation, without acid-induced denaturation.³⁸ Even the subsequent loss of solvent only causes relatively small changes. The final gas-phase structures still resemble the initial solution conformations, albeit with a slight compaction (Fig. 7C and Fig. S22C, S23C, S24C).^100,101,103 Thus, our findings bolster the view that solution-like structures can survive the transition from solution into the gas phase during native ESI,^{11,12,14–18}

Our simulations support the CRM/Rayleigh framework, which offers a simple way to predict native ESI charge states.^21,24,49,93 The model entails protein release into the gas phase upon droplet evaporation to dryness. Because the droplets stay close to the Rayleigh limit throughout their life cycle, the “residue” formed when the last solvent disappears has a net charge close to that of a protein-sized aqueous droplet. We demonstrate that this concept applies regardless of protein size and physicochemical parameters, for a pI range of at least 3.2 to 10.5. In addition, this model applies to native ESI in both positive and negative ion mode, provided that subtle stability differences of NH₄⁺-charged vs. Ac⁻-charged systems during the final stages of desolvation are taken into account.

Our findings contrast the proposal⁶⁷ that protein charging is governed by residue basicity. Support for our view that basicity is not critical for native ESI charging comes from experiments on proteins that behave according the CRM/Rayleigh framework, despite not possessing any titratable side chains (such that there is no PT between protein and solvent during ESI).⁹³ Native ESI charging for such proteins is mediated by NH₄⁺ adducts, similar to the pepsin behavior seen here. Pertinent insights also come from MD data on aqueous ESI droplets containing Na⁺ as the only type of charge carrier. [M + zNa]^z+ ions formed from such droplets have charge states matching those of experimentally observed [M + zH]^z+ ions.⁶⁵ This consistent charging behavior is despite the fact that protonation and Na⁺ binding are very different types of events, the latter not being governed by basicity.

Earlier investigations highlighted possible alternatives to the CRM/Rayleigh framework, demonstrating that small proteins may undergo IEM ejection.²⁵ However, protein IEM is viable only for compact conformers that carry a large charge (e.g. ubiquitin¹²⁺ at pH 3).²⁵ For the conditions studied here, droplets were buffered at pH 4.76 ± 1, thereby precluding the accumulation of a large solution charge on ubiquitin. It is therefore unsurprising that the current simulations did not show protein IEM. Instead, we observed CRM behavior regardless of protein and droplet size. The persistence of compact structures under native ESI conditions also precludes ion formation via the chain ejection mechanism (CEM), a pathway pursued by unfolded proteins.⁶⁵

The streamlined dMPMD algorithm used here allows ESI droplet simulations in an unprecedented size regime. By tracking the behavior of these large systems through the entire sequence of events (fission, evaporation, gas-phase ion formation), we were able to capture a much larger segment of the overall ESI process than previous modeling studies. In the future, we plan to extend this range even further, all the way to the ESI capillary and the formation of ESI droplets at the Taylor cone. Initial progress toward this goal has been achieved,¹⁰⁷ albeit not for macromolecular analytes such as proteins.

Author contributions

K. H. and L. K. conceived this project, developed computer code, and wrote the manuscript. K. H. performed all simulations and experiments.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting the findings of this study are available within the article and its supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6cp00481d.

The code is available on request from the authors.

Acknowledgements

Funding was provided by the Natural Sciences and Engineering Research Council of Canada (RGPIN-2024-05244).

References

1. J. B. Fenn, Angew. Chem., Int. Ed., 2003, 42, 3871–3894.
2. C. Kirschbaum, J. L. Bennett and C. V. Robinson, Anal. Chem., 2025, 97, 19331–19339.
3. Y. Yang, C. D. Niu, C. E. Bobst and I. A. Kaltashov, Anal. Chem., 2021, 93, 3337–3342.
4. H. Osterholz, A. Stevens, M. L. Abramsson, D. Lama, K. Brackmann, A. Rising, A. Elofsson, E. G. Marklund, S. Deindl, A. Leppert and M. Landreh, JACS Au, 2025, 5, 281–290.
5. A. G. Anders, E. D. Tidwell, V. V. Gadkari, M. Koutmos and B. T. Ruotolo, J. Am. Chem. Soc., 2024, 146, 4412–4420.
6. D. T. Bui, E. N. Kitova, P. I. Kitov, L. Han, L. K. Mahal and J. S. Klassen, J. Am. Chem. Soc., 2024, 146, 28809–28821.
7. S. Tamara, M. A. den Boer and A. J. R. Heck, Chem. Rev., 2022, 122, 7269–7326.
8. D. Cubrilovic, K. Barylyuk, D. Hofmann, M. J. Walczak, M. Graber, T. Berg, G. Wider and R. Zenobi, Chem. Sci., 2014, 5, 2794–2803.
9. J. Bellamy-Carter, M. Mohata, M. Falcicchio, J. Basran, Y. Higuchi, R. G. Doveston and A. C. Leney, Chem. Sci., 2021, 12, 10724–10731.
10. M. W. Zhou, C. Lantz, K. A. Brown, Y. Ge, L. Pasa-Tolic, J. A. Loo and F. Lermyte, Chem. Sci., 2020, 11, 20.
11. T. C. Cropley, F. C. Liu, M. Q. Chai, M. F. Bush and C. Bleiholder, J. Am. Chem. Soc., 2024, 146, 11115–11125.
12. J. L. P. Benesch and C. V. Robinson, Nature, 2009, 462, 576–577.
13. C. C. Moore, V. N. Staroverov and L. Konermann, J. Phys. Chem. B, 2023, 127, 4061–4071.
14. M. S. Westphall, K. W. Lee, A. Z. Salome, J. M. Lodge, T. Grant and J. J. Coon, Nat. Commun., 2022, 13, 6.
15. T. K. Esser, J. Böhning, P. Fremdling, M. T. Agasid, A. Costin, K. Fort, A. Konijnenberg, J. Gilbert, A. Bahm, A. Makarov, C. V. Robinson, J. L. P. Benesch, L. Baker, T. A. M. Bharat, J. Gault and S. Rauschenbach, PNAS Nexus, 2022, 1, 1–13.
16. D. E. Clemmer, D. H. Russell and E. R. Williams, Accounts Chem. Res., 2017, 50, 556–560.
17. T. M. Allison, E. Reading, I. Liko, A. J. Baldwin, A. Laganowsky and C. V. Robinson, Nat. Commun., 2015, 6, 8551.
18. S. M. Dixit, D. A. Polasky and B. T. Ruotolo, Curr. Opin. Chem. Biol., 2018, 42, 93–100.
19. R. V. Williams, C. Huang, K. W. Moremen, I. J. Amster and J. H. Prestegard, Sci. Rep., 2022, 12, 14769.
20. K. Breuker and F. W. McLafferty, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 18145–18152.
21. P. Kebarle and U. H. Verkerk, Mass Spectrom. Rev., 2009, 28, 898–917.
22. R. L. Grimm and J. L. Beauchamp, J. Phys. Chem. A, 2010, 114, 1411–1419.
23. V. S. Avadhani, C. C. Harper, Z. M. Miller and E. R. Williams, J. Am. Chem. Soc., 2025, 147, 18853–18863.
24. J. Fernandez de la Mora, Anal. Chim. Acta, 2000, 406, 93–104.
25. E. Aliyari and L. Konermann, Anal. Chem., 2020, 92, 10807–10814.
26. J. V. Iribarne and B. A. Thomson, J. Chem. Phys, 1976, 64, 2287–2294.
27. I. G. Loscertales and J. Fernandez de la Mora, J. Chem. Phys, 1995, 103, 5041–5060.
28. K. Hanifi, P. M. Scrosati and L. Konermann, J. Phys. Chem B, 2024, 128, 5973–5986.
29. S. Nguyen and J. B. Fenn, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 1111–1117.
30. N. Khristenko, F. Rosu, E. Largy, J. Haustant, C. Mesmin and V. Gabelica, J. Am. Chem. Soc., 2023, 145, 498–506.
31. R. R. Ogorzalek Loo, R. Lakshmanan and J. A. Loo, J. Am. Soc. Mass Spectrom., 2014, 25, 1675–1693.
32. J. Li, C. Santambrogio, S. Brocca, G. Rossetti, P. Carloni and R. Grandori, Mass Spectrom. Rev., 2016, 35, 111–122.
33. E. Christofi and P. Barran, Chem. Rev., 2023, 123, 2902–2949.
34. L. Konermann, E. Ahadi, A. D. Rodriguez and S. Vahidi, Anal. Chem., 2013, 85, 2–9.
35. S. Consta and J. K. Chung, J. Phys. Chem. B, 2011, 115, 10447–10455.
36. Z. Xia and E. R. Williams, J. Am. Soc. Mass Spectrom., 2018, 29, 194–202.
37. H. J. Sterling, C. A. Cassou, A. C. Susa and E. R. Williams, Anal. Chem., 2012, 84, 3795–3801.
38. A. Kharlamova, B. M. Prentice, T.-Y. Huang and S. A. McLuckey, Anal. Chem., 2010, 82, 7422–7429.
39. W. Wang, E. N. Kitova and J. S. Klassen, Anal. Chem., 2003, 75, 4945–4955.
40. M. M. Gadzuk-Shea, E. E. Hubbard, T. A. Gozzo and M. F. Bush, J. Am. Soc. Mass Spectrom., 2023, 34, 1675–1684.
41. G. Joseph, B. Binny and A. R. Venter, Anal. Chem., 2025, 97, 6034–6040.
42. Z. B. Han, N. Omata, T. Matsuda, S. Hishida, S. Takiguchi, R. Komori, R. Suzuki and L. C. Chen, Chem. Sci., 2023, 14, 4506–4515.
43. L. Konermann, Z. Y. Liu, Y. Haidar, M. J. Willans and N. A. Bainbridge, Anal. Chem., 2023, 95, 13957–13966.
44. T. E. Walker, A. Laganowsky and D. H. Russell, Anal. Chem., 2022, 94, 10824–10831.
45. X. Zhuang, A. F. M. Gavriilidou and R. Zenobi, J. Am. Soc. Mass Spectrom., 2017, 28, 341–346.
46. B. T. V. Davis, A. Velyvis and S. Vahidi, Anal. Chem., 2023, 95, 17525–17532.
47. J. R. Catai, J. S. Toraño, G. J. de Jong and G. W. Somsen, Analyst, 2007, 132, 75–81.
48. M. L. Kramer, H. D. Kratzin, B. Schmidt, A. Römer, O. Windl, S. Liemann, S. Hornemann and H. Kretzschmar, J. Biol. Chem., 2001, 276, 16711–16719.
49. I. A. Kaltashov and A. Mohimen, Anal. Chem., 2005, 77, 5370–5379.
50. R. J. Williams and C. M. Lyman, J. Am. Chem. Soc., 1932, 54, 1911–1912.
51. B. Naiman, J. Chem. Ed., 1948, 25, 454.
52. G. J. Van Berkel and V. Kertesz, Anal. Chem., 2007, 79, 5511–5520.
53. J. Y. Li, W. P. Lyu, G. Rossetti, A. Konijnenberg, A. Natalello, E. Ippoliti, M. Orozco, F. Sobott, R. Grandori and P. Carloni, J. Phys. Chem. Lett., 2017, 8, 1105–1112.
54. V. Znamenskiy, I. Marginean and A. Vertes, J. Phys. Chem. A, 2003, 107, 7406–7412.
55. C. Caleman and D. van der Spoel, Phys. Chem. Chem. Phys., 2007, 9, 5105–5111.
56. D. Kim, N. Wagner, K. Wooding, D. E. Clemmer and D. H. Russell, J. Am. Chem. Soc., 2017, 139, 2981–2988.
57. S. G. Kondalaji, M. Khakinejad and S. J. Valentine, J. Am. Soc. Mass Spectrom., 2018, 29, 1665–1677.
58. E. I. Calixte, O. T. Liyanage, H. J. Kim, E. D. Ziperman, A. J. Pearson and E. S. Gallagher, J. Phys. Chem. B, 2020, 124, 479–486.
59. H. Higashi, T. Tokumi, C. J. Hogan, H. Suda, T. Seto and Y. Otani, Phys. Chem. Chem. Phys., 2015, 17, 15746–15755.
60. R. Beveridge, L. G. Migas, R. K. Das, R. V. Pappu, R. W. Kriwacki and P. E. Barran, J. Am. Chem. Soc., 2019, 141, 4908–4918.
61. M. Porrini, F. Rosu, C. Rabin, L. Darre, H. Gomez, M. Orozco and V. Gabelica, ACS Central Sci., 2017, 3, 454–461.
62. B. Wang and D. P. Tieleman, Commun. Chem., 2023, 6, 21.
63. M. J. Luan, Z. H. Hou, B. Zhang, L. Ma, S. Yuan, Y. Liu and G. Huang, Anal. Chem., 2023, 95, 8798–8806.
64. A. Patriksson, E. Marklund and D. van der Spoel, Biochemistry, 2007, 46, 933–945.
65. L. Konermann, H. Metwally, Q. Duez and I. Peters, Analyst, 2019, 144, 6157–6171.
66. L. Konermann and S. Kim, J. Chem. Theo. Comp., 2022, 18, 3781–3794.
67. M. S. Cordes and E. S. Gallagher, J. Am. Chem. Soc., 2025, 147, 15066–15076.
68. S. K. Fegan and M. Thachuk, J. Chem. Theory Comput., 2013, 9, 2531–2539.
69. S. Nash and R. W. Vachet, J. Am. Chem. Soc., 2022, 144, 22128–22139.
70. J. Huang and A. D. MacKerell, J. Comput. Chem., 2013, 34, 2135–2145.
71. L. Konermann, H. Metwally, R. G. McAllister and V. Popa, Methods, 2018, 144, 104–112.
72. J. Abramson, J. Adler, J. Dunger, R. Evans, T. Green, A. Pritzel, O. Ronneberger, L. Willmore, A. J. Ballard, J. Bambrick, S. W. Bodenstein, D. A. Evans, C. C. Hung, M. O'Neill, D. Reiman, K. Tunyasuvunakool, Z. Wu, A. Zemgulyte, E. Arvaniti, C. Beattie, O. Bertolli, A. Bridgland, A. Cherepanov, M. Congreve, A. I. Cowen-Rivers, A. Cowie, M. Figurnov, F. B. Fuchs, H. Gladman, R. Jain, Y. A. Khan, C. M. R. Low, K. Perlin, A. Potapenko, P. Savy, S. Singh, A. Stecula, A. Thillaisundaram, C. Tong, S. Yakneen, E. D. Zhong, M. Zielinski, A. Zídek, V. Bapst, P. Kohli, M. Jaderberg, D. Hassabis and J. M. Jumper, Nature, 2024, 630, 24.
73. T. R. Covey, B. A. Thomson and B. B. Schneider, Mass Spectrom. Rev., 2009, 28, 870–897.
74. C. Vega and E. de Miguel, J. Chem. Phys, 2007, 126, 154707.
75. U. A. Mirza and B. T. Chait, Int. J. Mass. Spectrom. Ion Proc., 1997, 162, 173–181.
76. T. P. Silverstein, Front. Mol. Biosci., 2021, 8, 764099.
77. P. Mitchell, Nature, 1967, 214, 400.
78. M. Bakhtiari and L. Konermann, J. Phys. Chem. B, 2019, 123, 1784–1796.
79. G. A. Kaminski, R. A. Friesner, J. Tirado-Rives and W. L. Jorgensen, J. Phys. Chem. B, 2001, 105, 6474–6487.
80. A. Moser, K. Range and D. M. York, J. Phys. Chem. B, 2010, 114, 13911–13921.
81. T. Wyttenbach and M. T. Bowers, J. Phys. Chem. B, 2011, 115, 12266–12275.
82. S. J. Allen, A. M. Schwartz and M. F. Bush, Anal. Chem., 2013, 85, 12055–12061.
83. A. J. R. Heck and R. H. H. Van den Heuvel, Mass Spectrom. Rev., 2004, 23, 368–389.
84. C. Lantz, R. L. Rider, S.-T. Kuo, Z. Xi, E. Burningham, S. D. Yun, A. Laganowsky and D. H. Russell, J. Am. Chem. Soc., 2025, 147, 35616–35626.
85. B. T. Ruotolo, S.-J. Hyung, P. M. Robinson, K. Giles, R. H. Bateman and C. V. Robinson, Angew. Chem., Int. Ed., 2007, 46, 8001–8004.
86. P. D. Schnier, D. S. Gross and E. R. Williams, J. Am. Soc. Mass Spectrom., 1995, 6, 1086–1097.
87. A. G. Harrison, J. Am. Soc. Mass Spectrom., 2012, 23, 116–123.
88. S. Zhou and K. D. Cook, J. Am. Soc. Mass Spectrom., 2001, 12, 206–214.
89. A. P. Null, A. I. Nepomuceno and D. C. Muddiman, Anal. Chem., 2003, 75, 1331–1339.
90. Q. Duez, H. Metwally and L. Konermann, Anal. Chem., 2018, 90, 9912–9920.
91. S. Consta, J. Phys. Chem. B, 2022, 126, 8350–8357.
92. M. Dole, L. L. Mack, R. L. Hines, R. C. Mobley, L. D. Ferguson and M. B. Alice, J. Chem. Phys, 1968, 49, 2240–2249.
93. M. L. Abramsson, C. Sahin, J. T. S. Hopper, R. M. M. Branca, J. Danielsson, M. M. Xu, S. A. Chandler, N. Osterlund, L. L. Ilag, A. Leppert, J. Costeira-Paulo, L. S. Lang, K. Teilum, A. Laganowsky, J. L. P. Benesch, M. Oliveberg, C. V. Robinson, E. G. Marklund, T. M. Allison, J. R. Winther and M. Landreh, JACS Au, 2021, 1, 2385–2393.
94. G. Wang and R. B. Cole, Org. Mass Spectrom., 1994, 29, 419–427.
95. A. R. Dongré, J. L. Jones, Á. Somogyi and V. H. Wysocki, J. Am. Chem. Soc., 1996, 118, 8365–8374.
96. D. A. Polasky, S. M. Dixit, M. F. Keating, V. V. Gadkari, P. C. Andrews and B. T. Ruotolo, J. Am. Chem. Soc., 2020, 142, 6750–6760.
97. R. R. Abzalimov and I. A. Kaltashov, Anal. Chem., 2010, 82, 942–950.
98. N. Felitsyn, E. N. Kitova and J. S. Klassen, Anal. Chem., 2001, 73, 4647–4661.
99. Y. Mao, J. Woenckhaus, J. Kolafa, M. A. Ratner and M. F. Jarrold, J. Am. Chem. Soc., 1999, 121, 2712–2721.
100. Z. Hall, A. Politis, M. F. Bush, L. J. Smith and C. V. Robinson, J. Am. Chem. Soc., 2012, 134, 3429–3438.
101. A. D. Rolland, L. S. Biberic and J. S. Prell, J. Am. Soc. Mass Spectrom., 2022, 33, 369–381.
102. S.-H. Chen and D. H. Russell, J. Am. Soc. Mass Spectrom., 2015, 26, 1433–1443.
103. L. J. Persson, C. Sahin, M. Landreh and E. G. Marklund, Anal. Chem., 2024, 96, 15023–15030.
104. A. Patriksson, C. M. Adams, F. Kjeldsen, R. A. Zubarev and D. van der Spoel, J. Phys. Chem. B, 2007, 111, 13147–13150.
105. J. G. Bonner, Y. A. Lyon, C. Nellessen and R. R. Julian, J. Am. Chem. Soc., 2017, 139, 10286–10293.
106. K. Breuker, S. Brüschweiler and M. Tollinger, Angew. Chem., Int. Ed., 2011, 50, 873–877.
107. M. Dolatkhah Ouch Bolagh and L. Konermann, J. Phys. Chem. B, 2025, 129, 4096–4109.

---