Position-specific kinetic isotope effects for nitrous oxide: a new expansion of the Rayleigh model

Abstract

Nitrous oxide (N₂O) is a potent greenhouse gas and the most significant anthropogenic ozone-depleting substance currently being emitted. A major source of anthropogenic N₂O emissions is the microbial conversion of fixed nitrogen species from fertilizers in agricultural soils. Thus, understanding the enzymatic mechanisms by which microbes produce N₂O has environmental significance. Measurement of the ¹⁵N $/$ ¹⁴N isotope ratios of N₂O produced by purified enzymes or axenic microbial cultures is a promising technique for studying N₂O biosynthesis. Typically, N₂O-producing enzymes combine nitrogen atoms from two identical substrate molecules (NO or NH₂OH). Position-specific isotope analysis of the central (N^α) and outer (N^β) nitrogen atoms in N₂O enables the determination of the individual kinetic isotope effects (KIEs) for N^α and N^β, providing mechanistic insight into the incorporation of each nitrogen atom. Previously, position-specific KIEs (and fractionation factors) were quantified using the Rayleigh distillation equation, i.e., via linear regression of ${\mathit{\delta }}^{\mathrm{15}}$N^α or ${\mathit{\delta }}^{\mathrm{15}}$N^β against [$-f\ln f/(\mathrm{1}-f)$], where $f$ is the fraction of substrate remaining in a closed system. This approach, however, is inaccurate for N^α and N^β because it does not account for fractionation at N^α affecting the isotopic composition of substrate available for incorporation into the $\mathit{\beta }$ position (and vice versa). Therefore, we developed a new expansion of the Rayleigh model that includes specific terms for fractionation at the individual N₂O nitrogen atoms. By applying this Expanded Rayleigh model to a variety of simulated N₂O synthesis reactions with different combinations of normal, inverse, and/or no KIEs at N^α and N^β, we demonstrate that our new model is both accurate and robust. We also applied this new model to two previously published datasets describing N₂O production from NH₂OH oxidation in a methanotroph culture (Methylosinus trichosporium) and N₂O production from NO by a purified Histoplasma capsulatum (fungal) P450 NOR, demonstrating that the Expanded Rayleigh model is a useful tool in calculating position-specific fractionation for N₂O synthesis.

1 Introduction

Nitrous oxide (N₂O) is a greenhouse gas with a global warming potential approximately 300 times greater than that of CO₂ (100-year time horizon) (Stocker et al., 2013). In addition, N₂O is the primary source of stratospheric NO_x and is therefore a significant ozone-depleting substance (Ravishankara et al., 2009). Levels of atmospheric N₂O are increasing by an average of $\sim$ 0.29 % annually (Lan et al., 2022), primarily due to rising anthropogenic N₂O emissions, which comprise approximately 40 % of annual global N₂O emissions (Tian et al., 2020). In the US, at least 70 % of annual anthropogenic N₂O emissions result from agricultural practices, such as the application of nitrogen-containing fertilizer (e.g., ammonium and nitrate) to the field. Soil microbes can metabolize these nitrogen-containing compounds through processes such as nitrification and denitrification, producing N₂O either as a side product (nitrification), intermediate (complete denitrification), or end product (incomplete denitrification) (Kuypers et al., 2018). Understanding the underlying mechanisms of N₂O production for each metabolic pathway is critical for tracing and ultimately mitigating environmental N₂O emissions.

Many different types of enzymes catalyze N₂O formation in microbes. Examples of N₂O-producing enzyme families include cytochrome P460 (Caranto et al., 2016) and hydroxylamine oxidoreductase (Yamazaki et al., 2014) in nitrifiers, “bacterial” nitric oxide reductases (cNOR, qNOR, etc.) found in bacterial and archaeal denitrifiers, and fungal nitric oxide reductase (P450 NOR) from fungal denitrifiers (Yamazaki et al., 2014; Yang et al., 2014; Lehnert et al., 2018) as well as detoxifying enzymes such as flavodiiron nitric oxide reductases (Romão et al., 2016). Distinguishing among different sources of N₂O in the environment and identifying the relative contribution of each source are challenging and require a combination of analytical methods. One useful technique for N₂O source apportionment is stable isotope analysis. This method takes advantage of the fact that many physical and chemical processes, such as enzyme-catalyzed reactions, generally have slightly different rate constants for different isotopes of the same element (e.g., ¹⁴N and ¹⁵N). In an isotopically sensitive process, the isotope that reacts faster will be more prevalent (enriched) in the product relative to the substrate. Different processes generally have distinct patterns of isotopic enrichment and thus have unique isotopic signatures.

Stable isotope analysis is particularly useful for studying N₂O synthesis because the central (N^α) and outer (N^β) nitrogen atoms in N₂O are non-equivalent and non-exchangeable, and the isotope ratios of these two nitrogen atoms can be determined individually (Toyoda and Yoshida, 1999). These position-specific isotope ratios can aid in identifying which enzyme produces or which enzymes produce N₂O in microbial cultures under varying growth conditions (Sutka et al., 2003, 2004, 2006, 2008; Toyoda et al., 2005; Frame and Casciotti, 2010; Yamazaki et al., 2014; Haslun et al., 2018). Additionally, because isotopic enrichment is typically dependent on the rate-limiting step(s) of the enzyme-catalyzed reaction, analyzing the isotopic signatures of individual enzymes can provide valuable insights into each enzyme's catalytic mechanism.

Isotopic enrichment for a specific process can be quantified by determining a parameter known as the fractionation factor (${\mathit{\alpha }}_{\mathrm{p}/\mathrm{s}}$). For reactions with first-order kinetics, ${\mathit{\alpha }}_{\mathrm{p}/\mathrm{s}}$ is equal to the rate constant of the heavy isotope ($k_{\mathrm{H}}$) divided by the rate constant for the light isotope ($k_{\mathrm{L}}$). Because ${\mathit{\alpha }}_{\mathrm{p}/\mathrm{s}}$ values are typically close to 1, ${\mathit{\alpha }}_{\mathrm{p}/\mathrm{s}}$ is often converted to its per mil equivalent, the enrichment factor (${\mathit{\epsilon }}_{\mathrm{p}/\mathrm{s}}$) (Mariotti et al., 1981).

1$${\mathit{\epsilon }}_{\mathrm{p}/\mathrm{s}}=\left({\mathit{\alpha }}_{\mathrm{p}/\mathrm{s}}-\mathrm{1}\right)\cdot \mathrm{1000}$$ Notably, ${\mathit{\alpha }}_{\mathrm{p}/\mathrm{s}}$ is also equal to the inverse of the kinetic isotope effect (KIE $=k_{\mathrm{L}}$/$k_{\mathrm{H}}$). The KIE and fractionation factor values for a specific reaction describe the isotopic preference of that reaction. Most commonly, the light isotope reacts faster than the heavy isotope (i.e., KIE $\mathit{>}\mathrm{1}$), which is designated as a normal KIE. Conversely, when the heavy isotope reacts faster (KIE $\mathit{<}\mathrm{1}$), the reaction is said to have an inverse KIE. If the heavy and light isotopes react at the same rate (KIE $=$ 1), the reaction has no kinetic isotope effect.

Mariotti et al. (1981) derived an approximation of the Rayleigh distillation equation (hereafter referred to as the (standard) Rayleigh equation) that is frequently used to determine enrichment factors in closed systems (Eq. 2) (Mariotti et al., 1981). 2$${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{p}}={\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{s}\mathrm{0}}-{\mathit{\epsilon }}_{\mathrm{p}/\mathrm{s}}\cdot {\displaystyle \frac{f\ln \left(f\right)}{\mathrm{1}-f}}$$ Here, ${\mathit{\delta }}^{\mathrm{15}}$N^p is the $\mathit{\delta }$ value of accumulated product (p) at time $t$, ${\mathit{\delta }}^{\mathrm{15}}$N^s0 is the $\mathit{\delta }$ value of the substrate at $t=$ 0, and $f$ is the fraction of substrate remaining at time $t$ ($f$ $=$ N^s $/$ N^s0). (The $\mathit{\delta }$ values are directly related to isotope ratios; see the “Isotope nomenclature” section below.) Notably, the Rayleigh equation is based on the assumption that the element of interest is transferred from the substrate to a single position (or equivalent positions) in the product.

The Rayleigh equation (Eq. 2) holds true for ${\mathit{\delta }}^{\mathrm{15}}$N^bulk (the average of ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β) because the overall reaction can be thought of as the unidirectional conversion of two identical substrate molecules (e.g., two NO or two NH₂OH molecules) to one product molecule with two nitrogen atoms (2S$\to$P₂). In the absence of side reactions, this means that ${\mathit{\delta }}^{\mathrm{15}}$N^bulk can be calculated at any point during the reaction (i.e., for any value of $f$) as long as the $\mathit{\delta }$ value of the initial substrate (${\mathit{\delta }}^{\mathrm{15}}$N^s0) and the enrichment factor (${\mathit{\epsilon }}_{\text{N-bulk}}$) are known. In practice, the ${\mathit{\epsilon }}_{\text{N-bulk}}$ value for a specific enzymatic reaction is found by plotting ${\mathit{\delta }}^{\mathrm{15}}$N^bulk versus [$-f\ln f/(\mathrm{1}-f)$]. Because it is assumed that changes in ${\mathit{\delta }}^{\mathrm{15}}$N^bulk are due solely to the bulk kinetic isotope effect of the reaction in question, the slope of this plot (determined by linear regression) is the bulk enrichment factor (${\mathit{\epsilon }}_{\text{N-bulk}}$), and the $y$ intercept corresponds to ${\mathit{\delta }}^{\mathrm{15}}$N^s0.

Unfortunately, the relationship between ${\mathit{\delta }}^{\mathrm{15}}$N^α, ${\mathit{\delta }}^{\mathrm{15}}$N^β, and their corresponding enrichment factors is more complex because the two N atoms are non-exchangeable and both steps that incorporate NO (or NH₂OH) are potentially isotopically sensitive. In other words, fractionation at the $\mathit{\alpha }$ position alters the isotopic ratio of the remaining substrate (${\mathit{\delta }}^{\mathrm{15}}{N}^s$) that is available for incorporation into the $\mathit{\beta }$ position, and vice versa. Thus, simply determining the change in ${\mathit{\delta }}^{\mathrm{15}}$N^α (or ${\mathit{\delta }}^{\mathrm{15}}$N^β) as a function of [$-f\ln f/(\mathrm{1}-f)$] over the course of the reaction (i.e., finding the slope of a Rayleigh plot via linear regression of ${\mathit{\delta }}^{\mathrm{15}}$N^α or ${\mathit{\delta }}^{\mathrm{15}}$N^β against [$-f\ln f/(\mathrm{1}-f)]$) is not an accurate method for calculating the enrichment factor for N^α (or N^β).

The complexity of calculating enrichment factors for a branched N₂O biosynthesis reaction can be demonstrated mathematically by substituting $\mathrm{1}/\mathrm{2}$ (${\mathit{\delta }}^{\mathrm{15}}$N${}^{\mathit{\alpha }}+{\mathit{\delta }}^{\mathrm{15}}$N^β) for ${\mathit{\delta }}^{\mathrm{15}}$N^p and $\mathrm{1}/\mathrm{2}$ (${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}+{\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$) for ${\mathit{\epsilon }}_{\mathrm{p}/\mathrm{s}}$ in the Rayleigh equation (Eq. 2), yielding Eq. (3). 3$${\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\alpha }}+{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\beta }}\right)={\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{s}\mathrm{0}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}+{\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}\right){\displaystyle \frac{f\ln \left(f\right)}{\mathrm{1}-f}}$$ Here, ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$ stand for the isotopic enrichment factors for the $\mathit{\alpha }$ and $\mathit{\beta }$ nitrogen atoms, respectively. Equation (3) can then be rearranged to solve for ${\mathit{\delta }}^{\mathrm{15}}$N^α (Eq. 4). 4$$\begin{array}{rl}{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\alpha }}& =\mathrm{2}{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{s}\mathrm{0}}+{\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}\left(-{\displaystyle \frac{f\ln \left(f\right)}{\mathrm{1}-f}}\right)\\ & +{\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}\left(-{\displaystyle \frac{f\ln \left(f\right)}{\mathrm{1}-f}}\right)-{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\beta }}\end{array}$$ The underlying difficulty in using Eq. (4) is that for a given value of $f$, ${\mathit{\delta }}^{\mathrm{15}}$N^α depends not only on ${\mathit{\delta }}^{\mathrm{15}}$N^s0 and ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}$, but also on ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$ and the current ${\mathit{\delta }}^{\mathrm{15}}$N^β value. Put another way, we have one equation and two unknowns (${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$).

To address this issue, we developed an expanded version of the Rayleigh model that includes specific terms for the individual nitrogen atoms in N₂O. To test the robustness of this new Expanded Rayleigh model, we generated simulated isotopic data for each of five different potential catalytic mechanisms for N₂O formation. For each N₂O synthesis scenario, the size and skewness of the simulated measurement error were varied, and 1000 simulated datasets were generated for each combination of error level and skewness type. We then applied the standard Rayleigh model and Expanded Rayleigh model to each dataset to determine ${\mathit{\epsilon }}_{\text{N-bulk}}$, ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}$, and ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$ and compared the precision and accuracy of the two models. Finally, we compared the results of applying the standard Rayleigh model and Expanded Rayleigh model to previously published isotopic data.

2 Methods

2.1 Isotope nomenclature

The isotopic ratio, $R$, is defined as the abundance of the heavy isotope (¹⁵N) divided by the abundance of the light isotope (¹⁴N).

5$$R={\displaystyle \frac{\mathrm{Heavy}\mathrm{isotope}}{\mathrm{Light}\mathrm{isotope}}}=\left(\mathrm{e}.\mathrm{g}.,{\displaystyle \frac{{}^{\mathrm{15}}\mathrm{N}}{{}^{\mathrm{14}}\mathrm{N}}}\right)$$ In the case of nitrogen, ¹⁴N is naturally much more abundant than ¹⁵N (the atom percent of ¹⁵N in atmospheric N₂ (air) is 0.3663 $\pm$ 0.0004 %, so $R$ $=$ 0.0036765) (Junk and Svec, 1958; Mariotti et al., 1981; Skrzypek and Dunn, 2020). The $R$ value of a sample is usually measured relative to that of a standard ($R_{\mathrm{sample}}/R_{\mathrm{standard}}$). For nitrogen, the standard is atmospheric N₂. Because there is naturally little variation in ¹⁵N, $R_{\mathrm{sample}}/R_{\mathrm{standard}}$ values are typically close to 1, and isotopic ratios are usually converted to the $\mathit{\delta }$ notation, which is commonly expressed in the per mil form (‰) by multiplying the unitless form of $\mathit{\delta }$ by 1000: 6$$\mathit{\delta }=\left(\left({\displaystyle \frac{R_{\mathrm{sample}}}{R_{\mathrm{standard}}}}\right)-\mathrm{1}\right)\cdot \mathrm{1000},$$ where $R_{\mathrm{sample}}$ is the isotopic ratio of the sample, and $R_{\mathrm{standard}}$ is the isotopic ratio of the analytical standard (Mariotti et al., 1981). The per mil form of $\mathit{\delta }$ is used throughout this paper.

N^bulk is defined as the total number of “equivalents” of nitrogen atoms in N₂O, or $\mathrm{2}\cdot$ (mol of N₂O). N^α and N^β are defined as the number of moles of N at the $\mathit{\alpha }$ (central) and $\mathit{\beta }$ (outer) positions of N₂O, respectively. 7$$\mathrm{0.5}\cdot {\mathrm{N}}^{\mathrm{bulk}}={\mathrm{N}}^{\mathit{\alpha }}={\mathrm{N}}^{\mathit{\beta }}=\mathrm{mol}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}$$ Thus, ¹⁵N^bulk is defined as the total number of moles of ¹⁵N at both the $\mathit{\alpha }$ and $\mathit{\beta }$ positions. 8$${}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}{=}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\alpha }}{+}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\beta }}$$ Similarly, ¹⁴N^bulk is defined as the total number of moles of ¹⁴N at the $\mathit{\alpha }$ and $\mathit{\beta }$ positions of N₂O. 9$${}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{bulk}}{=}^{\mathrm{14}}{\mathrm{N}}^{\mathit{\alpha }}{+}^{\mathrm{14}}{\mathrm{N}}^{\mathit{\beta }}$$

2.2 $\mathit{\alpha }$, $\mathit{\epsilon }$, and KIE

As discussed above, the isotopic fractionation factor (${\mathit{\alpha }}_{\mathrm{p}/\mathrm{s}}$) is equal to the inverse of the kinetic isotope effect (KIE) for first-order reactions and for higher-order reactions when the rare (heavy) isotope is dilute (i.e., ¹⁵N or ¹³C at natural abundance) (Eq. 10) (Bigeleisen and Wolfsberg, 1958). KIE is defined as the rate constant of the light isotope, $k_{\mathrm{L}}$, divided by the rate constant of the heavy isotope, $k_{\mathrm{H}}$ (Eq. 11). $$\begin{array}{}\text{10}& {\displaystyle }\mathrm{KIE}={\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{p}/\mathrm{s}}}}\text{11}& {\displaystyle }\mathrm{KIE}={\displaystyle \frac{k_{\mathrm{L}}}{k_{\mathrm{H}}}}\end{array}$$ Formally, ${\mathit{\alpha }}_{\mathrm{p}/\mathrm{s}}$ is defined as the isotopic ratio of instantaneously produced product ($R_{\mathrm{pi}}$) divided by the isotopic ratio of residual substrate at the same instant ($R_{\mathrm{s}}$) (Eq. 12) (Mariotti et al., 1981).

12$${\mathit{\alpha }}_{\mathrm{p}/\mathrm{s}}={\displaystyle \frac{\frac{{}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{pi}}}{{}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{pi}}}}{\frac{{}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{s}}}{{}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{s}}}}}={\displaystyle \frac{R_{\mathrm{pi}}}{R_{\mathrm{s}}}}$$ As noted above, the fractionation factor can also be expressed in per mil (‰) terms by using Eq. (1) to convert $\mathit{\alpha }$ to the per mil enrichment factor, $\mathit{\epsilon }$.

2.3 Assumptions

We used the Rayleigh equation (Eq. 2) to determine ${\mathit{\alpha }}_{\text{N-bulk}}$, the fractionation factor for N^bulk. The Rayleigh model depends on five key assumptions, as outlined below.

2.3.1 Assumption 1

The reaction is irreversible and occurs in a closed system in the absence of side reactions (Mariotti et al., 1981).

This means that the final ${\mathit{\delta }}^{\mathrm{15}}$N^bulk value is equal to the initial ${\mathit{\delta }}^{\mathrm{15}}$N^s value (${\mathit{\delta }}^{\mathrm{15}}$N^s0).

13$${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk},\mathrm{final}}={\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{s}\mathrm{0}}$$

2.3.2 Assumption 2

The bulk fractionation factor (${\mathit{\alpha }}_{\mathrm{p}/\mathrm{s}}$ or ${\mathit{\alpha }}_{\text{N-bulk}}$) is constant.

This means that the isotopic ratio of bulk N produced in an infinitely short time period ($R_{\mathrm{bulk},i}$) divided by the isotopic ratio of the remaining substrate ($R_{\mathrm{s}}$) is constant (Eq. 14) (Mariotti et al., 1981).

14$${\mathit{\alpha }}_{\text{N-bulk}}={\displaystyle \frac{R_{\mathrm{bulk},i}}{R_{\mathrm{s}}}}$$ This assumption should be valid as long as the rate-limiting step of the chemical transformation in question remains the same.

2.3.3 Assumption 3

The light isotope (¹⁴N) is much more abundant than the heavy isotope (¹⁵N), which is true for all experiments performed with ¹⁵N at or near natural abundance.

This assumption was employed by Mariotti and colleagues to approximate total N abundance (¹⁴N ${+}^{\mathrm{15}}$N) as ¹⁴N in their definition of $f$ (Mariotti et al., 1981). For N₂O isotopocules, the assumption that ¹⁴N $\gg$ ¹⁵N also means that variations in the apportionment of ¹⁵N between the $\mathit{\alpha }$ and $\mathit{\beta }$ positions will not significantly affect the abundance of ¹⁴N at the $\mathit{\alpha }$ or $\mathit{\beta }$ position. Thus, at or near the natural abundance of ¹⁵N, both ¹⁴N^α and ¹⁴N^β may be approximated as $\mathrm{0.5}{\cdot }^{\mathrm{14}}$N^bulk. This approximation underlies the assumption that ${\mathit{\delta }}^{\mathrm{15}}$N^bulk is approximately equal to the average of ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β (Eq. 15).

15$${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot \left({\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\alpha }}+{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\beta }}\right)$$ Equation (15) is routinely used in isotopic studies of N₂O to calculate ${\mathit{\delta }}^{\mathrm{15}}$N^β values using experimentally determined ${\mathit{\delta }}^{\mathrm{15}}$N^bulk and ${\mathit{\delta }}^{\mathrm{15}}$N^α values (Toyoda and Yoshida, 1999).

2.3.4 Assumption 4

The values of ${\mathit{\delta }}_{s\mathrm{0}}$ and ${\mathit{\delta }}_{\mathrm{s}}$ are relatively close to 0 ‰.

Mariotti's approximation of the Rayleigh equation uses a simplification that is only valid if ${\mathit{\delta }}_{\mathrm{s}}/\mathrm{1000}$ and ${\mathit{\delta }}_{s\mathrm{0}}/\mathrm{1000}$ are small relative to 1 (Mariotti et al., 1981). In practice, for an N₂O synthesis reaction where ${\mathit{\epsilon }}_{\text{N-bulk}}$ is $-$20 ‰, a ${\mathit{\delta }}^{\mathrm{15}}$N^s0 value of $\pm$ 100 ‰ will introduce an error of $\sim$ 2 ‰ when ${\mathit{\epsilon }}_{\text{N-bulk}}$ is determined using this approximation (Table S1).

2.3.5 Assumption 5

The absolute value of $\mathit{\epsilon }$ (e.g., ${\mathit{\epsilon }}_{\text{N-bulk}}$) is relatively small.

The magnitude of ${\mathit{\delta }}^{\mathrm{15}}$N^s values depends on both ${\mathit{\delta }}^{\mathrm{15}}$N^s0 and ${\mathit{\epsilon }}_{\text{N-bulk}}$. Thus, similar to the previous assumption, this assumption was introduced to facilitate the approximation that ${\mathit{\delta }}^{\mathrm{15}}$N${}^{\mathrm{s}}/\mathrm{1000}$ is small relative to 1 (Mariotti et al., 1981). As with Assumption 4, the error introduced by this approximation increases as the absolute value of $\mathit{\epsilon }$ increases. For a reaction where ${\mathit{\delta }}^{\mathrm{15}}$N^s0 is equal to 0 ‰ and ${\mathit{\epsilon }}_{\text{N-bulk}}$ is equal to $\pm$ 50 ‰, this approximation introduces an error of $\sim$ 1 ‰ to $\mathit{\epsilon }$ values (${\mathit{\epsilon }}_{\text{N-bulk}}$, ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}$, and ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$), while a reaction where ${\mathit{\epsilon }}_{\text{N-bulk}}$ is equal to $\pm$ 100 ‰ introduces an error of 5 ‰–6 ‰ (Table S2).

2.3.6 Assumption 6

For the purposes of modeling N₂O synthesis, we will use all of the well-established assumptions listed above. Additionally, to account for fractionation at N^α and N^β individually, we introduce one more assumption: the individual fractionation factors for the $\mathit{\alpha }$ and $\mathit{\beta }$ nitrogen atoms in N₂O (${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$) are essentially constant.

We define ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$ in a manner analogous to that of ${\mathit{\alpha }}_{\text{N-bulk}}$ (Eq. 14). Thus, ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ is defined as the isotopic ratio of N^α from instantaneously produced N₂O ($R_{\mathit{\alpha },i}$) divided by $R_{\mathrm{s}}$ (Eq. 16).

16$${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}={\displaystyle \frac{\frac{{}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\alpha },i}}{{}^{\mathrm{14}}{\mathrm{N}}^{\mathit{\alpha },i}}}{\frac{{}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{s}}}{{}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{s}}}}}={\displaystyle \frac{R_{\mathit{\alpha },i}}{R_{\mathrm{s}}}}$$ Similarly, ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$ is defined as the instantaneous isotopic ratio of N^β ($R_{\mathit{\beta },i}$) divided by $R_{\mathrm{s}}$ (Eq. 17). 17$${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}={\displaystyle \frac{\frac{{}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\beta },i}}{{}^{\mathrm{14}}{\mathrm{N}}^{\mathit{\beta },i}}}{\frac{{}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{s}}}{{}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{s}}}}}={\displaystyle \frac{R_{\mathit{\beta },i}}{R_{\mathrm{s}}}}$$ Put another way, we assume a priori that the processes that govern preferential incorporation of one isotope over another remain unchanged at both the $\mathit{\alpha }$ and $\mathit{\beta }$ positions of N₂O during the course of the reaction. We expect this assumption to be valid throughout the reaction, as long as the rate-limiting steps of N₂O biosynthesis are unchanged.

2.4 The standard Rayleigh model and simulation of error-free ${\mathit{\delta }}^{\mathrm{15}}$N^bulk and ${\mathit{\delta }}^{\mathrm{15}}$N^s values

To compare the standard Rayleigh approach and our nonlinear expansion of the Rayleigh model, we simulated isotopic data for a number of different combinations of kinetic isotope effects for N^α and N^β (see Table 1 and Fig. 1). Each idealized, error-free dataset had three groups of equally spaced values of $f$ (0.7–0.3), similar to a typical experiment with three replicates (Tables S4–8). For each dataset, ${\mathit{\delta }}^{\mathrm{15}}$N^s0 was set at 0 ‰, and ${\mathit{\epsilon }}_{\text{N-bulk}}$ was set at either $-$20 ‰ (normal KIE, KIE $=$ 1.0204) or $+$20 ‰ (inverse KIE, KIE $=$ 0.9804). Equation (2), Mariotti's approximation of the Rayleigh model for accumulated product (Mariotti et al., 1981), was used to calculate ${\mathit{\delta }}^{\mathrm{15}}$N^bulk values of N₂O for each value of $f$ (${\mathit{\delta }}^{\mathrm{15}}$N${}^{\mathrm{bulk}}={\mathit{\delta }}^{\mathrm{15}}$N^p). Similarly, $\mathit{\delta }$ values for remaining substrate (${\mathit{\delta }}^{\mathrm{15}}$N^s) were calculated using Eq. (18), Mariotti's approximation of the Rayleigh model written in terms of substrate (Mariotti et al., 1981).

18$${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{s}}={\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{s}\mathrm{0}}+{\mathit{\epsilon }}_{\text{N-bulk}}\cdot \ln \left(f\right)$$ The ${\mathit{\delta }}^{\mathrm{15}}$N^bulk and ${\mathit{\delta }}^{\mathrm{15}}$N^s values simulated using the standard Rayleigh equation are shown in Fig. 1 in gray and orange, respectively.

Table 1

Parameters for idealized, error-free datasets describing N₂O synthesis for five simulated scenarios (Datasets 1–5) with varying combinations of ${\mathit{\epsilon }}_{\text{N-bulk}}$ and $\mathit{\rho }$.

Simulated	Subset	${\mathit{\epsilon }}_{\text{N-bulk}}$	${\mathit{\delta }}_{\mathrm{s}\mathrm{0}}$	${\mathit{\rho }}^{\mathrm{b}}$	${\mathit{\tau }}^{\mathrm{c}}$
dataset		(‰)a	(‰)
1	0.3 $\le$ $f$ $\le$ 0.7	$-$20	0	0.5102	0.49996
2	0.3 $\le$ $f$ $\le$ 0.7	$-$20	0	0.5050	0.49998
3	0.3 $\le$ $f$ $\le$ 0.7	$-$20	0	0.5200	0.49993
4	0.3 $\le$ $f$ $\le$ 0.7	20	0	0.5050	0.49998
5	0.3 $\le$ $f$ $\le$ 0.7	$-$20	0	0.4898	0.50004

^a Note that an $\mathit{\epsilon }$ value of $-$20 ‰ corresponds to a KIE of 1.0204 and an $\mathit{\epsilon }$ value of $+$20 ‰ corresponds to a KIE of 0.9804. ^b Parameter introduced in the Expanded Rayleigh model (Eq. 23): $\mathit{\rho }$ $=$ ¹⁵N${}^{\mathit{\alpha }}{/}^{\mathrm{15}}$N^bulk. ^c Parameter introduced in the Expanded Rayleigh model (Eq. 24): $\mathit{\tau }$ $=$ ¹⁴N${}^{\mathit{\alpha }}{/}^{\mathrm{14}}$N^bulk. $\mathit{\tau }$ was calculated for each value of $f$; the average value of $\mathit{\tau }$ is listed.

Figure 1

Simulated ${\mathit{\delta }}^{\mathrm{15}}$N values for N₂O synthesis in a closed system with different combinations of KIE ¹⁵N^α and KIE ¹⁵N^β. ${\mathit{\delta }}^{\mathrm{15}}$N values of remaining substrate (${\mathit{\delta }}^{\mathrm{15}}$N^s) and accumulated product (N₂O) (${\mathit{\delta }}^{\mathrm{15}}$N^bulk, ${\mathit{\delta }}^{\mathrm{15}}$N^α, and ${\mathit{\delta }}^{\mathrm{15}}$N^β) are plotted as a function of the fraction of substrate reduced ($\mathrm{1}-f$) with reactions progressing from left to right. Values of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk and ${\mathit{\delta }}^{\mathrm{15}}$N^s were calculated using Eq. (2) or Eq. (18) (Mariotti's approximation of the Rayleigh equation; Mariotti et al., 1981), where ${\mathit{\epsilon }}^{\mathrm{15}}$N^bulk is set to $-$20 ‰ (Datasets 1–3 and 5) or $+$20 ‰ (Dataset 4) and the initial ${\mathit{\delta }}^{\mathrm{15}}$N^s value (${\mathit{\delta }}^{\mathrm{15}}$N^s0) is set to 0 ‰. Values of ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β were calculated using Eqs. (27) and (28). For Datasets 1 and 5, either ${\mathit{\epsilon }}^{\mathrm{15}}$N^α or ${\mathit{\epsilon }}^{\mathrm{15}}$N^β was set to 0 ‰ (no isotope effect, KIE $=$ 1), and the corresponding value of $\mathit{\rho }$ was back-calculated (see the Supplement for details). For Datasets 2–4, values of ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β were calculated using pre-determined values of $\mathit{\rho }$ that were less than (Datasets 2 and 4) or greater than (Dataset 3) the value of $\mathit{\rho }$ from Dataset 1. See Table 3 for the KIE values for each simulation.

[Figure omitted. See PDF]

2.5 Nonlinear expansion of the Rayleigh model

Simulated ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β values were calculated based on the assumption that, like the fractionation factor for bulk N, the fractionation factors for the individual N^α and N^β atoms are constant (see Assumption 6 above). To relate ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$ to ${\mathit{\alpha }}_{\text{N-bulk}}$, we define the fraction of instantaneously converted ¹⁵N^bulk (¹⁵N^bulk,i) that is incorporated into the $\mathit{\alpha }$ position as $\mathit{\rho }$ (Eq. 19).

19$$\mathit{\rho }={\displaystyle \frac{{}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\alpha },i}}{{}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk},i}}}$$ Similarly, we define the fraction of instantaneously converted ¹⁴N (¹⁴N^bulk,i) that goes to the $\mathit{\alpha }$ position as $\mathit{\tau }$ (Eq. 20). 20$$\mathit{\tau }={\displaystyle \frac{{}^{\mathrm{14}}{\mathrm{N}}^{\mathit{\alpha },i}}{{}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{bulk},i}}}$$ In keeping with the premise that ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$ are constant, we assume a priori that $\mathit{\rho }$ and $\mathit{\tau }$ are constant. Thus, the equation for ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ may be written in terms of $\mathit{\rho }$, $\mathit{\tau }$, and ${\mathit{\alpha }}_{\text{N-bulk}}$. 21$${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}={\displaystyle \frac{\frac{\mathit{\rho }{\cdot }^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk},i}}{\mathit{\tau }{\cdot }^{\mathrm{14}}{\mathrm{N}}^{\mathrm{bulk},i}}}{\frac{{}^{\mathrm{15}}{N}^{\mathrm{s}}}{{}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{s}}}}}={\displaystyle \frac{\mathit{\rho }}{\mathit{\tau }}}\cdot {\mathit{\alpha }}_{\text{N-bulk}}$$ In the absence of side reactions, ¹⁵N^bulk,i must be equal to the sum of ¹⁵N^α,i and ¹⁵N^β,i. Similarly, ¹⁴N^bulk,i must be equal to the sum of ¹⁴N^α,i and ¹⁴N^β,i. Thus, the equation for ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$ may also be written in terms of $\mathit{\rho }$, $\mathit{\tau }$, and ${\mathit{\alpha }}_{\text{N-bulk}}$. 22$${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}={\displaystyle \frac{\frac{\left(\mathrm{1}-\mathit{\rho }\right){\cdot }^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk},i}}{\left(\mathrm{1}-\mathit{\tau }\right){\cdot }^{\mathrm{14}}{\mathrm{N}}^{\mathrm{bulk},i}}}{\frac{{}^{\mathrm{15}}{N}_{\mathrm{s}}}{{}^{\mathrm{14}}{\mathrm{N}}_{\mathrm{s}}}}}={\displaystyle \frac{\mathrm{1}-\mathit{\rho }}{\mathrm{1}-\mathit{\tau }}}\cdot {\mathit{\alpha }}_{\text{N-bulk}}$$ Equations (21) and (22) form the basis for the Expanded Rayleigh model.

To simulate ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β values and determine the corresponding fractionation factors, $\mathit{\rho }$ and $\mathit{\tau }$ must be related to accumulated isotope ratios, $R_{\mathit{\alpha }}$ and $R_{\mathit{\beta }}$. Although $\mathit{\rho }$ is formally defined in Eq. (19) as the fraction of instantaneously converted ¹⁵N apportioned to the $\mathit{\alpha }$ position, $\mathit{\rho }$ can also be written in terms of accumulated ¹⁵N values as long as $\mathit{\rho }$ is constant. 23$$\mathit{\rho }={\displaystyle \frac{{}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\alpha }}}{{}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}}}$$ Similarly, as long as $\mathit{\tau }$ is constant, $\mathit{\tau }$ may be written in terms of accumulated ¹⁴N values. 24$$\mathit{\tau }={\displaystyle \frac{{}^{\mathrm{14}}{\mathrm{N}}^{\mathit{\alpha }}}{{}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{bulk}}}}$$ Using the definition of $\mathit{\rho }$ from Eq. (23), $R_{\mathit{\alpha }}$ can be written in terms of $\mathit{\rho }$ by substituting ¹⁵N^α with $\mathit{\rho }{\cdot }^{\mathrm{15}}$N^bulk and replacing ¹⁴N^α with 0.5 $\cdot$ (N^bulk) - ¹⁵N^α (Eq. 25). 25$$R_{\mathit{\alpha }}={\displaystyle \frac{{}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\alpha }}}{{}^{\mathrm{14}}{\mathrm{N}}^{\mathit{\alpha }}}}={\displaystyle \frac{\mathit{\rho }{\cdot }^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}}{\mathrm{0.5}\cdot {\mathrm{N}}^{\mathrm{bulk}}-\mathit{\rho }{\cdot }^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}}}$$ The same approach can be used to define $R_{\mathit{\beta }}$ in terms of $\mathit{\rho }$ (Eq. 26). 26$$R_{\mathit{\beta }}={\displaystyle \frac{{}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\beta }}}{{}^{\mathrm{14}}{\mathrm{N}}^{\mathit{\beta }}}}={\displaystyle \frac{(\mathrm{1}-\mathit{\rho }){\cdot }^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}}{\mathrm{0.5}\cdot {\mathrm{N}}^{\mathrm{bulk}}-(\mathrm{1}-\mathit{\rho }){\cdot }^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}}}$$ The term 0.5 $\cdot$ N^bulk is equal to moles of N₂O (Eq. 7), and ¹⁵N^bulk can be calculated for any value of $f$ using N^bulk and ${\mathit{\delta }}^{\mathrm{15}}$N^bulk (see the Supplement). Equations (25) and (26) can be converted to $\mathit{\delta }$ notation as follows (see Eq. 6). $$\begin{array}{}\text{27}& {\displaystyle }{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\alpha }}=\left({\displaystyle \frac{\mathit{\rho }{\cdot }^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}}{(\mathrm{0.5}\cdot {\mathrm{N}}^{\mathrm{bulk}}-\mathit{\rho }{\cdot }^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}})\cdot R_{\mathrm{standard}}}}-\mathrm{1}\right)\cdot \mathrm{1000}\text{28}& {\displaystyle }{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\beta }}=\left({\displaystyle \frac{\left(\mathrm{1}-\mathit{\rho }\right){\cdot }^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}}{(\mathrm{0.5}\cdot {\mathrm{N}}^{\mathrm{bulk}}-(\mathrm{1}-\mathit{\rho }\left){\cdot }^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}\right)\cdot {\mathrm{R}}_{\mathrm{standard}}}}-\mathrm{1}\right)\cdot \mathrm{1000}\end{array}$$ Equations (27) and (28) can be used to determine ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β based on ${\mathit{\delta }}^{\mathrm{15}}$N^bulk as long as $\mathit{\rho }$ is constant and ${\mathit{\delta }}^{\mathrm{15}}$N^bulk can be approximated as the average of ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β (see Eq. 15 and Table S3).

2.6 Simulation of error-free ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β values

Idealized, error-free ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β values (shown in Fig. 1) were simulated for the values of $\mathit{\rho }$ listed in Table 1 by calculating ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β using Eqs. (27) and (28) (Tables S4–8). As outlined in Table 1, five different scenarios with different combinations of isotope effects (none, normal, or inverse) for the $\mathit{\alpha }$ and $\mathit{\beta }$ nitrogen atoms of N₂O were simulated. For Datasets 1 and 5, either ${\mathit{\epsilon }}^{\mathrm{15}}$N^α or ${\mathit{\epsilon }}^{\mathrm{15}}$N^β was set to 0 ‰ (no isotope effect, KIE $=$ 1), and the corresponding values of $\mathit{\rho }$ and $\mathit{\tau }$ were back-calculated (see the Supplement for details). The remaining values of $\mathit{\rho }$ were chosen to produce datasets with the different combinations of normal and/or inverse isotope effects (Datasets 2–4). The $\mathit{\rho }$ values tested (0.4898–0.5200) are all close enough to 0.5 such that the average of ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β is approximately equal to ${\mathit{\delta }}^{\mathrm{15}}$N^bulk (Table S3). In all cases, discrepancies between ${\mathit{\delta }}^{\mathrm{15}}$N^bulk values calculated using Eq. (2) and ${\mathit{\delta }}^{\mathrm{15}}$N^bulk values calculated by averaging ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β (Eq. 15) are well below the typical analytical error of 0.5 ‰–0.7 ‰ (Yang et al., 2014).

2.7 Determination of $\mathit{\rho }$ and $\mathit{\tau }$ for the Expanded Rayleigh model

To determine ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$ using Eqs. (21) and (22), $\mathit{\rho }$ and $\mathit{\tau }$ must be calculated using values that can be measured experimentally (see Fig. 2 for the general workflow for using the Expanded Rayleigh model to determine ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$). Because the abundance of ¹⁴N is so much greater than that of ¹⁵N, variations in ¹⁵N apportionment between N^α and N^β are negligible when determining the fraction of ¹⁴N^bulk at the $\mathit{\alpha }$ position (i.e., $\mathit{\tau }$). For enzyme-catalyzed N₂O synthesis with substrate at or near the natural abundance of ¹⁵N, we have found that $\mathit{\tau }$ is always 0.5000 (see the analysis of data from Sutka et al., 2006, and Yang et al., 2014, below). For the purposes of testing our new Expanded Rayleigh model, we calculated more precise values of $\mathit{\tau }$. To do this, $\mathit{\tau }$ was calculated for each value of $f$ using Eq. (24), and the average $\mathit{\tau }$ value for each dataset was used to calculate fractionation factors.

To determine $\mathit{\rho }$ using values that can be measured experimentally, Eq. (15) is rewritten so that $\mathit{\rho }$ is the only unknown parameter. This can be accomplished by replacing ${\mathit{\delta }}^{\mathrm{15}}$N^α with the equivalent term from Eq. (27), producing Eq. (29), a nonlinear equation we have named nonlinear model 1.

29$$\begin{array}{rl}{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot [\left({\displaystyle \frac{\mathit{\rho }{\cdot }^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}}{\left(\mathrm{0.5}\cdot {\mathrm{N}}^{\mathrm{bulk}}-\mathit{\rho }{\cdot }^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}\right)\cdot R_{\mathrm{standard}}}}-\mathrm{1}\right)\\ & \cdot \mathrm{1000}+{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\beta }}]\end{array}$$ Alternatively, ${\mathit{\delta }}^{\mathrm{15}}$N^β could be replaced with a term that includes $\mathit{\rho }$ (Eq. 28), producing Eq. (30) (nonlinear model 2). 30$$\begin{array}{rl}{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\cdot [{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\alpha }}\\ & +\left({\displaystyle \frac{\left(\mathrm{1}-\mathit{\rho }\right){\cdot }^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}}{\left(\mathrm{0.5}\cdot {\mathrm{N}}^{\mathrm{bulk}}-(\mathrm{1}-\mathit{\rho }){\cdot }^{\mathrm{15}}{\mathrm{N}}^{\mathrm{bulk}}\right)\cdot R_{\mathrm{standard}}}}-\mathrm{1}\right)\\ & \cdot \mathrm{1000}]\end{array}$$ Both nonlinear models were used to find the optimal value of $\mathit{\rho }$ via nonlinear least-squares regression for each of our simulated datasets.

2.8 Generation of datasets with simulated error

To mimic experimental error, 1000 simulated datasets were generated from each error-free dataset by adding random error to N^s, ${\mathit{\delta }}^{\mathrm{15}}$N^bulk, and ${\mathit{\delta }}^{\mathrm{15}}$N^α using a modified version of the procedure outlined by Scott and colleagues (Scott et al., 2004). Error was propagated to $f$ by recalculating $f$ using the N^s values with error added (f $=$ N^s / N^s0). Similarly, error was propagated to ${\mathit{\delta }}^{\mathrm{15}}$N^β by recalculating ${\mathit{\delta }}^{\mathrm{15}}$N^β using ${\mathit{\delta }}^{\mathrm{15}}$N^bulk and ${\mathit{\delta }}^{\mathrm{15}}$N^α values with error added (Eq. 15). For each simulation, fractionation factors and KIE values were calculated using the standard Rayleigh model and the Expanded Rayleigh model (versions 1 and 2), and the KIE values from each set of 1000 simulations were averaged. Three “levels” of error (low, medium, and high) and three types of skewness (none, left-skewed, and right-skewed; see below) were used for a total of nine conditions (Table 2; Figs. 3 and S1–4). For simulations with low error, the standard deviation of each randomly generated error term was set to the level of error expected in a typical experiment (Table 2). The standard deviation of N^s values was set to 1.5 % of N^s0 (the maximum value of N^s), which corresponds to the estimated level of error in gas volume measurements (see the Supplement). The standard deviations of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk and ${\mathit{\delta }}^{\mathrm{15}}$N^α were set to 0.5 ‰ and 0.7 ‰, respectively, which correspond to standard deviations typical for isotopic analysis of N₂O via isotope ratio mass spectrometry (IRMS) (Yang et al., 2014). To further test the robustness of the linear and nonlinear models, medium- and high-level error terms were generated by multiplying the low-level error terms for N^s, ${\mathit{\delta }}^{\mathrm{15}}$N^bulk, and ${\mathit{\delta }}^{\mathrm{15}}$N^α by the appropriate factor (Table 2). The multiplication factors for medium- and high-level error were chosen so that the average $R^{\mathrm{2}}$ values produced by linear regression of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk against [$-f\ln f/(\mathrm{1}-f)$] were approximately 0.7 and 0.4, respectively. Similar $R^{\mathrm{2}}$ values have been observed for standard Rayleigh linear regression plots for N₂O biosynthesis (e.g., Yang et al., 2014; Figs. S2–4).

Figure 2

General workflow for using the Expanded Rayleigh model to determine fractionation factors ($\mathit{\alpha }$) and kinetic isotope effects (KIEs) for the individual N atoms (N^α and N^β) in N₂O. To determine the fractionation factors for the $\mathit{\alpha }$ and $\mathit{\beta }$ atoms in N₂O (${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$), three separate values must be calculated. (1) The fractionation factor for bulk N (${\mathit{\alpha }}_{\text{N-bulk}}$) is determined by finding the slope (${\mathit{\epsilon }}_{\text{N-bulk}}$) of a standard Rayleigh plot where $y$ $=$ [$-f\ln f/(\mathrm{1}-f)$] and $x={\mathit{\delta }}^{\mathrm{15}}$N^bulk and converting ${\mathit{\epsilon }}_{\text{N-bulk}}$ to ${\mathit{\alpha }}_{\text{N-bulk}}$ as has been previously described (Mariotti et al., 1981). (2) Nonlinear modeling is used to determine the optimal value of $\mathit{\rho }$ (the fraction of ¹⁵N^bulk at the $\mathit{\alpha }$ position). For Expanded Rayleigh model 1 (shown here), ${\mathit{\delta }}^{\mathrm{15}}$N^α is replaced in the top equation with an equivalent expression that includes $\mathit{\rho }$ (bottom equation) to generate the model used for nonlinear least-squares regression (Eq. 29). For Expanded Rayleigh model 2, ${\mathit{\delta }}^{\mathrm{15}}$N^β is substituted, producing Eq. (30) (not shown). (3) The value of $\mathit{\tau }$ (the fraction of ¹⁴N^bulk at the $\mathit{\alpha }$ position) is determined by dividing ¹⁴N^α by ¹⁴N^bulk for each step of the reaction and averaging the results. (4) Finally, ${\mathit{\alpha }}_{\text{N-bulk}}$, $\mathit{\rho }$, and $\mathit{\tau }$ are used to calculate ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$. (5–6) These individual fractionation factors can be converted to the corresponding KIE values (Eq. 10) or $\mathit{\epsilon }$ values (Eq. 1).

[Figure omitted. See PDF]

Figure 3

Example simulations derived from Dataset 1 with varying levels of error and types of skewness. Each panel shows a single dataset (representative of 1000 simulated datasets) consisting of three replicates with five time points each. All simulations were derived from Dataset 1 (${\mathit{\delta }}^{\mathrm{15}}$N${}^{\mathrm{s}\mathrm{0}}=$ 0 ‰, ${\mathit{\epsilon }}_{\text{N-bulk}}=$ $-$20 ‰ (KIE ¹⁵N${}^{\mathrm{bulk}}=$ 1.0204), KIE ¹⁵N${}^{\mathit{\alpha }}=$ 1.0000, and KIE ¹⁵N${}^{\mathit{\beta }}=$ 1.0417). Each level of error and type of skewness is described in Table 2. Example graphs derived from Datasets 2–4 are shown in Figs. S1–4.

[Figure omitted. See PDF]

For each level of error, three sets of 1000 simulations with different types of skewness (no skew, left-skewed, or right-skewed) were performed to mimic types of skewness that might be observed in experimental data. For example, if N₂O leaks from a sample bottle prior to analysis, the quantity of N₂O produced will be underestimated, and the corresponding values of N^s and $f$ will be overestimated. Additionally, since light isotopes have a higher rate of both diffusion and effusion, values of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk, ${\mathit{\delta }}^{\mathrm{15}}$N^α, and ${\mathit{\delta }}^{\mathrm{15}}$N^β will be overestimated. This scenario is represented by simulations where the random error values added to N^s, ${\mathit{\delta }}^{\mathrm{15}}$N^bulk, and ${\mathit{\delta }}^{\mathrm{15}}$N^α were derived from left-skewed distributions. Skewness classifications were then assigned based on the skewness of the residuals for the standard Rayleigh model where the dependent variable ($y$) is ${\mathit{\delta }}^{\mathrm{15}}$N^bulk (Table 2). Simulations with residuals between $-$0.5 and 0.5 (inclusive) were designated “no skew,” meaning the distribution of these residuals did not differ significantly from the normal distribution. Simulations where standard Rayleigh residuals had a skewness value less than $-$0.5 were designated “left-skewed” (i.e., left-tailed distribution), and simulations where skewness was greater than 0.5 were designated “right-skewed” (i.e., right-tailed distribution).

Table 2

Standard deviation and skewness of random error added to simulated N^s, ${\mathit{\delta }}^{\mathrm{15}}$N^bulk, and ${\mathit{\delta }}^{\mathrm{15}}$N^α values.

Level of errora	Ns (nmol)	${\mathit{\delta }}^{\mathrm{15}}$Nbulk (‰)	${\mathit{\delta }}^{\mathrm{15}}$Nα (‰)	Skewnessb
Low	0.015 $\cdot$ Ns0c	0.5	0.7	None	$|$Skewness$|$ $\le$ 0.5
				Left	Skewness $\mathit{<}-\mathrm{0.5}$
				Right	Skewness $\mathit{>}\mathrm{0.5}$
Medium	2 $\cdot$ (0.015 $\cdot$ Ns0)c	2.5 $\cdot$ (0.5)	2.5 $\cdot$ (0.7)	None	$|$Skewness$|$ $\le$ 0.5
				Left	Skewness $\mathit{<}-\mathrm{0.5}$
				Right	Skewness $\mathit{>}\mathrm{0.5}$
High	3 $\cdot$ (0.015 $\cdot$ Ns0)c	5 $\cdot$ (0.5)	5 $\cdot$ (0.7)	None	$|$Skewness$|$ $\le$ 0.5
				Left	Skewness $\mathit{<}-\mathrm{0.5}$
				Right	Skewness $\mathit{>}\mathrm{0.5}$

^a “Level of error” refers to the standard deviations of N^s (nmol of remaining substrate), ${\mathit{\delta }}^{\mathrm{15}}$N^bulk, and ${\mathit{\delta }}^{\mathrm{15}}$N^a. ^b Skewness of the residuals for the standard Rayleigh model (linear regression of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk against [$-f\ln f/(\mathrm{1}-f)$]). ^c N^s0 refers to nanomoles (nmol) of substrate at time 0 (i.e., the maximum amount of substrate). N^s0 was set to 10 000 nmol in all simulations.

Table 3

KIE values for N^bulk, N^α, and N^β calculated directly from simulation input values (Actual) or from the standard Rayleigh model or Expanded Rayleigh model applied to error-free simulated datasets.

Simulated dataset		Actuala	Model
			Standard	Expanded	Expanded
			Rayleighb	Rayleigh 1c	Rayleigh 2c
1	KIE 15Nbulk	1.0204	1.0204	1.0204	1.0204
(No KIE 15Nα,	KIE 15Nα	1.0000	1.0208	1.0000	1.0000
Normal KIE 15Nβ)	KIE 15Nβ	1.0417	1.0200	1.0417	1.0417
2	KIE 15Nbulk	1.0204	1.0204	1.0204	1.0204
(Normal KIE 15Nα,	KIE 15Nα	1.0103	1.0206	1.0103	1.0103
Normal KIE 15Nβ)	KIE 15Nβ	1.0308	1.0202	1.0308	1.0308
3	KIE 15Nbulk	1.0204	1.0204	1.0204	1.0204
(Inverse KIE 15Nα,	KIE 15Nα	0.9810	1.0212	0.9810	0.9810
Normal KIE 15Nβ)	KIE 15Nβ	1.0631	1.0196	1.0631	1.0631
4	KIE 15Nbulk	0.9804	0.9804	0.9804	0.9804
(Inverse KIE 15Nα,	KIE 15Nα	0.9706	0.9802	0.9707	0.9706
Inverse KIE 15Nβ)	KIE 15Nβ	0.9903	0.9806	0.9903	0.9903
5	KIE 15Nbulk	1.0204	1.0204	1.0204	1.0204
(Normal KIE 15Nα,	KIE 15Nα	1.0417	1.0200	1.0417	1.0417
No KIE 15Nβ)	KIE 15Nβ	1.0000	1.0208	1.0000	1.0000

^a Actual values were calculated directly from the input values used to generate the dataset. ^b KIE values were calculated from ${\mathit{\epsilon }}_{\text{N-bulk}}$, ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}$, or ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$ values obtained via linear regression of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk, ${\mathit{\delta }}^{\mathrm{15}}$N^α, or ${\mathit{\delta }}^{\mathrm{15}}$N^β against [$-f\ln f/(\mathrm{1}-f)$] (see Eq. 2). ^c For the Expanded Rayleigh model, ${\mathit{\alpha }}_{\text{N-bulk}}$ was determined with the standard Rayleigh approach, $\mathit{\rho }$ was determined via nonlinear regression (nonlinear model 1 or 2, Eq. 29 or Eq. 30), and $\mathit{\tau }$ was determined by averaging ¹⁴N${}^{\mathit{\alpha }}{/}^{\mathrm{14}}$N^bulk for every step of the reaction.

To assess the accuracy of each model under different circumstances, average KIE values and standard deviations were calculated for each set of 1000 simulations with error. The absolute relative difference of each average KIE was calculated using Eq. (31): 31$$\text{absolute relative difference}=\mathrm{|}{\displaystyle \frac{\text{estimate}-\text{actual}}{\text{actual}}}\mathrm{|},$$ where the estimate is the average KIE and actual values were calculated directly from the input values for simulations without error (see Table 3). Lower absolute relative difference values indicate higher accuracy. Goodness of fit was assessed for the (linear) standard Rayleigh model and nonlinear models 1 and 2 by calculating the average root mean square error (RMSE) for each set of 1000 simulations. For each simulation, $p$ values for ${\mathit{\epsilon }}_{\text{N-bulk}}$ and $\mathit{\rho }$ (coefficients extracted from the (linear) standard Rayleigh model and nonlinear model 1 or 2, respectively) were determined using the one-sample $t$ test (two-sided) (Baty et al., 2015; Ritz and Streibig, 2008; Kalpiæ et al., 2011). The null hypothesis used to calculate these $p$ values is that there is no kinetic isotope effect. For the linear model, this means that the null hypothesis is ${\mathit{\epsilon }}_{\text{N-bulk}}=$ 0 (KIE ¹⁵N${}^{\mathrm{bulk}}=$ 1). For the nonlinear models, the null hypothesis is that $\mathit{\rho }=$ 0.5 (KIE ¹⁵N${}^{\mathit{\alpha }}=$ 1 and KIE ¹⁵N${}^{\mathit{\beta }}=$ 1).

2.9 Analysis of previously published experimental data

The standard and Expanded Rayleigh models were applied to published data on N₂O produced from NH₂OH by Methylosinus trichosporium (ATCC 49243) (Sutka et al., 2006). We calculated values of $f$, ${\mathit{\delta }}^{\mathrm{15}}$N^α, and ${\mathit{\delta }}^{\mathrm{15}}$N^β using the N₂O concentrations, ${\mathit{\delta }}^{\mathrm{15}}$N^bulk values, and site preference (SP) values from replicate B measured by Sutka et al. (see Table S14 for original data and our calculated values and Fig. 4 for a graphical representation of these data). SP is defined as the difference between ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β.

32$$\mathrm{SP}={\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\alpha }}-{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}^{\mathit{\beta }}$$ The data for replicate B were chosen because these data had a range of $f$ values large enough for ${\mathit{\epsilon }}_{\text{N-bulk}}$ to be determined via the standard Rayleigh method (linear regression of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk against [$-f\ln f/(\mathrm{1}-f)$]). See the Supplement for details.

Figure 4

Experimentally measured ${\mathit{\delta }}^{\mathrm{15}}$N values for N₂O synthesis from NH₂OH in a pure culture of Methylosinus trichosporium (data modified from Sutka et al., 2006). Data previously reported for M. trichosporium (Methylocystis sp.) replicate B (Sutka et al., 2006). ${\mathit{\delta }}^{\mathrm{15}}$N values of accumulated product (N₂O) (${\mathit{\delta }}^{\mathrm{15}}$N^bulk, ${\mathit{\delta }}^{\mathrm{15}}$N^α, and ${\mathit{\delta }}^{\mathrm{15}}$N^β) are plotted as a function of the fraction of substrate consumed ($\mathrm{1}-f$) with the reaction progressing from left to right. We back-calculated values of $f$, ${\mathit{\delta }}^{\mathrm{15}}$N^α, and ${\mathit{\delta }}^{\mathrm{15}}$N^β using the reported N₂O concentrations, ${\mathit{\delta }}^{\mathrm{15}}$N^bulk values, and SP values. The initial measured ${\mathit{\delta }}^{\mathrm{15}}$N value of NH₂OH (${\mathit{\delta }}^{\mathrm{15}}$N^s0), $-$2.3 ‰, is shown at $f$ $=$ 1.

[Figure omitted. See PDF]

We also applied the standard and Expanded Rayleigh models to previously published isotopic data on N₂O production from purified Histoplasma capsulatum (fungal) P450 NOR (Fig. 5 and Table S16) (Yang et al., 2014). Because this dataset covered a wide range of $f$ values ($f$ $=$ 0.42-0.87) and ${\mathit{\delta }}^{\mathrm{15}}$N^bulk varied linearly with [$-f\ln f/(\mathrm{1}-f)$], we used the standard Rayleigh model (Eq. 2) to determine ${\mathit{\epsilon }}_{\text{N-bulk}}$ via linear regression of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk against [$-f\ln f/(\mathrm{1}-f)$]. ${\mathit{\epsilon }}_{\text{N-bulk}}$ was determined for each of three independent biological replicates (13 observations per replicate), and the results were averaged. Similarly, to determine standard Rayleigh values of ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$, we performed linear regression of ${\mathit{\delta }}^{\mathrm{15}}$N^α or ${\mathit{\delta }}^{\mathrm{15}}$N^β against [$-f\ln f/(\mathrm{1}-f)$] for each independent replicate and averaged the results.

Figure 5

Experimentally measured product ${\mathit{\delta }}^{\mathrm{15}}$N values and calculated ${\mathit{\delta }}^{\mathrm{15}}$N^s0 values for N₂O synthesis from NO by purified H. capsulatum P450 NOR. Adapted with permission from Yang et al. (2014). Copyright 2014 American Chemical Society. Data previously reported for H. capsulatum (fungal) P450 NOR replicates A–C (Yang et al., 2014). ${\mathit{\delta }}^{\mathrm{15}}$N values of accumulated product (N₂O) (${\mathit{\delta }}^{\mathrm{15}}$N^bulk, ${\mathit{\delta }}^{\mathrm{15}}$N^α, and ${\mathit{\delta }}^{\mathrm{15}}$N^β) are plotted as a function of the fraction of substrate reduced ($\mathrm{1}-f$) with the reaction progressing from left to right. We back-calculated ${\mathit{\delta }}^{\mathrm{15}}$N^s0 for each replicate by finding the intercept of the corresponding standard Rayleigh plot (Eq. 2) where $x$ $=$ [$-f\ln f/(\mathrm{1}-f)$] and $y={\mathit{\delta }}^{\mathrm{15}}$N^bulk.

[Figure omitted. See PDF]

To apply the Expanded Rayleigh model to the isotopic data on N₂O production by P450 NOR, we used a modified procedure because, in this case, $\mathit{\rho }$ and $\mathit{\tau }$ vary linearly with [$-f\ln f/(\mathrm{1}-f)$]. Using nonlinear regression to determine $\mathit{\rho }$ yields a value of $\mathit{\rho }$ in between the more extreme values of $\mathit{\rho }$ observed at the start and end of the reaction and does not represent the data from the overall reaction very well. Therefore, instead of using nonlinear model 1 (Eq. 29) or model 2 (Eq. 30) to determine $\mathit{\rho }$ using data from the entire reaction, we calculated $\mathit{\rho }$ for each individual time point (i.e., each individual value of $f$) using Eq. (23). Similarly, $\mathit{\tau }$ was calculated for each individual observation using Eq. (24). These individual $\mathit{\rho }$ and $\mathit{\tau }$ values were then used to calculate ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ (Eq. 21) and ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$ (Eq. 22) values for each value of $f$, which were then converted to KIE values using Eq. (10). This approach yielded KIE values that increased (KIE ¹⁵N^α) or decreased (KIE ¹⁵N^β) as the reaction progressed. To estimate the range of position-specific KIEs produced by the Expanded Rayleigh model, individual KIE ¹⁵N^α and KIE ¹⁵N^β values representing the early part of the reaction ($f$ $=$ 0.8 $\pm$ 0.03) or the latter part of the reaction (late reaction, $f$ $=$ 0.5 $\pm$ 0.03) were averaged (Tables 6 and S17).

2.10 Modeling, statistical analysis, and figures

Modeling and statistical analyses of simulated and experimental $\mathit{\delta }$ values were performed with R statistical software (R Core Team, 2022), and figures were produced with ggplot2 (Wickham, 2016). To determine $\mathit{\rho }$ for the Expanded Rayleigh model, nonlinear least-squares regression was performed as previously described (Baty et al., 2015) using a starting $\mathit{\rho }$ value of 0.5.

For datasets with simulated error, random numbers representing simulated error were generated using the rsn function from the skew-normal distribution package (Azzalini, 2023), and skewness was calculated with the moments package (Komsta and Novomestky, 2022). In the moments package, skewness is defined as $(\mathrm{1}/n)\cdot \mathrm{\Sigma }\left(\right(x-\stackrel{\mathrm{‾}}{x})/s{)}^{\mathrm{3}}$, where $n$ is the sample size, $\stackrel{\mathrm{‾}}{x}$ is the sample mean, and $s$ is the sample standard deviation (Hippel, 2011).

For previously published experimental data, a linear model was used to determine if SP, $\mathit{\rho }$, or $\mathit{\tau }$ varied as a function of [$-f\ln f/(\mathrm{1}-f)$].

3 Results

3.1 Error-free simulations

To demonstrate that the standard Rayleigh model produces inaccurate results for the individual nitrogen atoms in N₂O, idealized, error-free datasets were simulated representing different scenarios with varying combinations of KIEs for N^α and N^β (Fig. 1 and Table 1). Because these data were simulated assuming that the fractions of ¹⁵N and ¹⁴N apportioned to each position remain constant (i.e., constant $\mathit{\rho }$ and $\mathit{\tau })$, the distance between ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β over the course of each reaction is constant, and the rates of change in ${\mathit{\delta }}^{\mathrm{15}}$N^α, ${\mathit{\delta }}^{\mathrm{15}}$N^β, and ${\mathit{\delta }}^{\mathrm{15}}$N^bulk with respect to reaction progress ($\mathrm{1}-f$) are essentially equal (Fig. 1). Thus, when the standard Rayleigh model (Eq. 2) is applied to each dataset, the slopes ($\mathit{\epsilon })$ of ${\mathit{\delta }}^{\mathrm{15}}$N^α, ${\mathit{\delta }}^{\mathrm{15}}$N^β, and ${\mathit{\delta }}^{\mathrm{15}}$N^bulk against [$-f\ln f/(\mathrm{1}-f)$] for each dataset are all approximately equal, as are the corresponding KIE ¹⁵N^bulk, KIE ¹⁵N^α, and KIE ¹⁵N^β values (Table 3). While the standard Rayleigh KIE ¹⁵N^bulk values match the actual KIE ¹⁵N^bulk values, the KIE ¹⁵N^α and KIE ¹⁵N^β values determined using the standard Rayleigh approach differ significantly from the actual KIEs calculated from simulation input values (Table 3). In each of the five simulated reactions (Datasets 1–5), the isotopic preference for N^α differs from that of N^β, which can be verified visually by noting that the ${\mathit{\delta }}^{\mathrm{15}}$N^α values are significantly different than the ${\mathit{\delta }}^{\mathrm{15}}$N^β values throughout each reaction (Fig. 1). However, in all five cases, the standard Rayleigh model produces KIE values for N^α and N^β that are approximately equal, highlighting the fact that the standard Rayleigh model inaccurately quantifies ¹⁵N apportionment between N^α and N^β. (If KIE ¹⁵N^α were equal to KIE ¹⁵N^β, the curves for ${\mathit{\delta }}^{\mathrm{15}}$N^α, ${\mathit{\delta }}^{\mathrm{15}}$N^β, and ${\mathit{\delta }}^{\mathrm{15}}$N^bulk shown in Fig. 1 would all be on top of each other.)

The inadequacies of the standard Rayleigh model are most evident when N^α and N^β have opposing kinetic isotope effects, i.e., scenarios where there is an inverse isotope effect at one position and a normal isotope effect at the other position (Dataset 3), or when only one position has a nonzero fractionation factor (Datasets 1 and 5). In these cases, the standard Rayleigh model fails to predict the correct type of isotope effect for one of the nitrogen atoms. For example, Dataset 1 represents a scenario where there is no isotopic preference at the $\mathit{\alpha }$ position, meaning instantaneous ${\mathit{\delta }}^{\mathrm{15}}$N^α (${\mathit{\delta }}^{\mathrm{15}}$N^αi) is always equal to ${\mathit{\delta }}^{\mathrm{15}}$N^s, which can be verified visually by noting that accumulated ${\mathit{\delta }}^{\mathrm{15}}$N${}^{\mathit{\alpha }}\approx {\mathit{\delta }}^{\mathrm{15}}$N^s at the start of the reaction (Fig. 1). Under these circumstances, KIE ¹⁵N^α is equal to 1.0000. However, the standard Rayleigh model predicts a much higher KIE ¹⁵N^α value (1.0208), incorrectly indicating that there is a reasonably strong, normal KIE at the $\mathit{\alpha }$ position. The standard Rayleigh KIE ¹⁵N^β value calculated for Dataset 1 is also too low (1.0200 instead of 1.0417), underestimating the preference for ¹⁴N at the $\mathit{\beta }$ position. Similarly misleading results are produced for Dataset 5, where there is no isotope effect for N^β, but the standard Rayleigh model predicts a normal isotope effect at this position. Along the same lines, the standard Rayleigh KIE values for Dataset 3 incorrectly indicate that both N^α and N^β have normal isotope effects when, in this case, N^α actually has an inverse isotope effect. Thus, even for a simulated dataset with no experimental error, the standard Rayleigh model produces results that are quantitatively – and sometimes qualitatively – incorrect.

The standard Rayleigh model produces inaccurate results for the individual N atoms in N₂O because the changes in ${\mathit{\delta }}^{\mathrm{15}}$N^α depend on both KIE ¹⁵N^α and KIE ¹⁵N^β, as do changes in ${\mathit{\delta }}^{\mathrm{15}}$N^β. Linear regression of standard Rayleigh plots for ${\mathit{\delta }}^{\mathrm{15}}$N^α or ${\mathit{\delta }}^{\mathrm{15}}$N^β fails to capture this nuance. To address this issue, we developed the Expanded Rayleigh model that introduces two new parameters, $\mathit{\rho }$ and $\mathit{\tau }$, to define how ¹⁵N^bulk and ¹⁴N^bulk, respectively, are apportioned between N^α and N^β. Fractionation factors for N^α and N^β are determined by combining ${\mathit{\alpha }}_{\text{N-bulk}}$, $\mathit{\rho }$, and $\mathit{\tau }$ (Eqs. 21 and 22). For experiments conducted at the natural abundance of ¹⁵N, the value of $\mathit{\tau }$, which is equal to ¹⁴N^α/¹⁴N^bulk, should always be quite close to 0.5, as discussed above. Indeed, each simulated reaction step has a mean value of $\mathit{\tau }$ equal to 0.500 for Datasets 1–5 (Tables S9–13), so the average value of ¹⁴N${}^{\mathit{\alpha }}{/}^{\mathrm{14}}$N^bulk was used to calculate ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$. To determine the value of $\mathit{\rho }$ that fits the entire dataset as accurately as possible, we performed nonlinear least-squares regression using the formula shown in Eq. (29) (nonlinear model 1) or Eq. (30) (nonlinear model 2). For error-free Datasets 1–5, nonlinear least-squares regression with either nonlinear model correctly predicted $\mathit{\rho }$. Combining these $\mathit{\rho }$ and $\mathit{\tau }$ values with ${\mathit{\alpha }}_{\text{N-bulk}}$ values derived from the standard Rayleigh equation yields KIE ¹⁵N^α and KIE ¹⁵N^β values that are identical to the expected values out to three or four decimal places (Table 3). Thus, unlike the standard Rayleigh approach, the Expanded Rayleigh model successfully recapitulates the kinetic isotope effects for N^α and N^β in an error-free dataset.

Table 4

Precision and accuracy of values calculated with Expanded Rayleigh model 1 using simulated datasets derived from Dataset 1 (no isotope effect for N^α, normal isotope effect for N^β). n/a – not applicable.

Model	Level of	Skewness	Standard	KIE 15Nbulka	KIE 15Nαa	KIE 15Nβa
	error		Rayleigh $R^{{\mathrm{2}}^{\mathrm{a},\mathrm{b}}}$
Actualc	None	No	n/a	1.0204 $\pm$ n/a	1 $\pm$ n/a	1.0417 $\pm$ n/a
Expanded	Low	No	0.95 $\pm$ 0.02	1.0202 $\pm$ 0.0013	0.9998 $\pm$ 0.0013	1.0414 $\pm$ 0.0014
Rayleigh 1d		Left	0.94 $\pm$ 0.02	1.0202 $\pm$ 0.0013	0.9998 $\pm$ 0.0013	1.0414 $\pm$ 0.0014
		Right	0.94 $\pm$ 0.02	1.0202 $\pm$ 0.0013	0.9998 $\pm$ 0.0013	1.0415 $\pm$ 0.0014
Expanded	Medium	No	0.75 $\pm$ 0.09	1.0196 $\pm$ 0.0031	0.9992 $\pm$ 0.0030	1.0408 $\pm$ 0.0032
Rayleigh 1d		Left	0.74 $\pm$ 0.1	1.0197 $\pm$ 0.0032	0.9994 $\pm$ 0.0032	1.0409 $\pm$ 0.0033
		Right	0.74 $\pm$ 0.1	1.0197 $\pm$ 0.0031	0.9993 $\pm$ 0.0031	1.0410 $\pm$ 0.0033
Expanded	High	No	0.44 $\pm$ 0.17	1.0186 $\pm$ 0.0060	0.9982 $\pm$ 0.0060	1.0398 $\pm$ 0.0061
Rayleigh 1d		Left	0.43 $\pm$ 0.17	1.0190 $\pm$ 0.0058	0.9987 $\pm$ 0.0058	1.0401 $\pm$ 0.0061
		Right	0.43 $\pm$ 0.17	1.0187 $\pm$ 0.0059	0.9983 $\pm$ 0.0059	1.0400 $\pm$ 0.0062

^a Average value $\pm$ standard deviation calculated from 1000 values generated from 1000 simulated datasets. ^b Average $R^{\mathrm{2}}$ value for linear regression of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk against [$-f\ln f/(\mathrm{1}-f)$] (using Eq. 2, the standard Rayleigh equation). ^c Actual values were calculated directly from the input values used to generate the dataset. ^d For Expanded Rayleigh model 1, ${\mathit{\alpha }}_{\text{N-bulk}}$ was determined with the standard Rayleigh approach, $\mathit{\rho }$ was determined via nonlinear regression (nonlinear model 1, Eq. 29), and $\mathit{\tau }$ was determined by averaging ¹⁴N${}^{\mathit{\alpha }}{/}^{\mathrm{14}}$N^bulk for every step of the reaction (Eq. 24). The results for Expanded Rayleigh model 2 (Table S9) are very similar.

3.2 Application of the Expanded Rayleigh model to simulations with error

To test the robustness of the Expanded Rayleigh model, we applied this model to simulated data with error at varying levels of size and skewness (Table 2; Figs. 3 and S1–4) and averaged the results from 1000 simulations for each error category (Tables 4 and S9–13). These results indicate that at a low level of error (the level of error expected for a typical experiment), the KIE values produced by the Expanded Rayleigh model are quite accurate (absolute difference relative to true values is 0.001, or 0.1 %, see Eq. 31), regardless of the skewness of the standard Rayleigh model residuals (Figs. 6 and S5–8). The standard deviations for each average KIE ¹⁵N^bulk, KIE ¹⁵N^α, or KIE ¹⁵N^β value (0.0013–0.0014) are also similar to standard error values from typical experiments ($\sim$ 0.0004–0.0021) (Tables 5 and 6) (Sutka et al., 2006; Yang et al., 2014). As the absolute value of the simulated error increases, the average KIE values become less accurate, and the corresponding standard deviations increase. Nonetheless, even at the highest level of error, the maximum absolute relative difference of average KIE values is 0.005 (0.5 % difference), and the range of values covered by average KIE $\pm$ standard deviation includes the actual KIE value. Thus, the Expanded Rayleigh model provides a robust method for determining KIE values even when the datasets have a high level of error or skewness.

Figure 6

Comparison of the accuracy and goodness of fit of the standard Rayleigh model and Expanded Rayleigh models 1 and 2. (a) Comparison of the accuracy of KIE ¹⁵N^α and KIE ¹⁵N^β values. For each model, the absolute relative difference for KIE ¹⁵N^α and KIE ¹⁵N^β values (average for 1000 simulated datasets derived from Dataset 1) is shown. Actual values for Dataset 1: ${\mathit{\delta }}^{\mathrm{15}}$N${}^{\mathrm{s}\mathrm{0}}=$ 0 ‰, ${\mathit{\epsilon }}_{\text{N-bulk}}=$ $-$20 ‰ (KIE ¹⁵N${}^{\mathrm{bulk}}=$ 1.0204), KIE ¹⁵N${}^{\mathit{\alpha }}=$ 1.0000, and KIE ¹⁵N${}^{\mathit{\beta }}=$ 1.0417. Absolute relative difference (Eq. 31) is the absolute value of the difference between the estimated value and actual value divided by the actual value ($|$(estimate – actual)/actual$|$). (b) Comparison of the average RMSE values for each set of 1000 simulated datasets derived from Dataset 1. Both the standard Rayleigh model and the Expanded Rayleigh model (1 and 2) use ${\mathit{\delta }}^{\mathrm{15}}$N^bulk as the dependent variable, so average RMSE values can be compared directly. Lower RMSE values indicate a better goodness of fit. At each error level (low, medium, and high), the absolute relative difference value (a) or average RMSE (b) is depicted with a symbol that represents skewness type as shown in the legend. Note that in most cases these symbols overlap.

[Figure omitted. See PDF]

In contrast to the Expanded Rayleigh model, the KIE ¹⁵N^α and KIE ¹⁵N^β values determined for simulations with error using the standard Rayleigh model alone are much less accurate (Figs. 6 and S5–8). For all simulations derived from Datasets 1–5 at all levels of error and types of skewness, the absolute relative differences of average standard Rayleigh KIE ¹⁵N^α and KIE ¹⁵N^β values ranged from 0.01–0.04 (1 %–4 % difference). These absolute relative differences are 2–40 times higher than the corresponding differences of Expanded Rayleigh model estimates (Figs. 6 and S5–8). Differences of this magnitude are large enough that the true KIE ¹⁵N^α and KIE ¹⁵N^β values typically do not fall within the range covered by average KIE $\pm$ standard deviation calculated using the standard Rayleigh model (data not shown). Thus, as discussed for the no-error simulations, standard Rayleigh KIE ¹⁵N^α and KIE ¹⁵N^β values do not accurately quantify isotopic fractionation at the $\mathit{\alpha }$ and $\mathit{\beta }$ positions of N₂O and can even lead to inaccurate designation of the type of KIE (i.e., no KIE, normal KIE, or inverse KIE) at either N^α or N^β. In short, for all levels of error and types of skewness tested, the standard Rayleigh model is only accurate for bulk N. To obtain KIE ¹⁵N^α and KIE ¹⁵N^β values with the same level of accuracy as KIE ¹⁵N^bulk, the KIE ¹⁵N^bulk value obtained via the standard Rayleigh model must be combined with $\mathit{\rho }$ and $\mathit{\tau }$ parameters from the Expanded Rayleigh model.

The error for KIE ¹⁵N^α and KIE ¹⁵N^β values is derived from three sources: error in ${\mathit{\alpha }}_{\text{N-bulk}}$, error in $\mathit{\rho }$, and error in $\mathit{\tau }$. In all of our simulations with added error, the most significant source of error is the error associated with ${\mathit{\alpha }}_{\text{N-bulk}}$, which is derived from linear regression of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk against [$-f\ln f/(\mathrm{1}-f)$]. As the level of error increases, the accuracy of ${\mathit{\epsilon }}_{\text{N-bulk}}$ (and thus the accuracy of ${\mathit{\alpha }}_{\text{N-bulk}}$ and KIE ¹⁵N^bulk) decreases more significantly than the accuracy of $\mathit{\rho }$ and $\mathit{\tau }$ estimates (Tables S9–13). The same trend is apparent for other assessments of model fit. The RMSE values for the nonlinear regression portion of the Expanded Rayleigh model (1 or 2) are always slightly lower than the RMSE values for the standard Rayleigh linear model (Figs. 6 and S5–8), indicating that the linear and nonlinear models produce similar “goodness of fit” values. Similarly, the $p$ values for ${\mathit{\epsilon }}_{\text{N-bulk}}$, the slope of the linear regression model, are always greater than the $p$ values for $\mathit{\rho }$, the coefficient extracted from the nonlinear portion of the Expanded Rayleigh model (in both cases, the null hypothesis is that KIE $=$ 1; see the Methods section). Out of each set of 1000 simulations at high-level error, a small number of simulated datasets (14 %–19 %) produce linear regression plots where confidence in the slope is low (i.e., ${\mathit{\epsilon }}_{\text{N-bulk}}$ $p$ value $\mathit{>}\mathrm{0.05}$) (data not shown). In contrast, all of the simulations derived from all five datasets with varying types of KIEs produce $\mathit{\rho }$ values that are likely to be significant ($p$ value $\mathit{<}\mathrm{0.001}$). Overall, application of the Expanded Rayleigh model to simulated datasets with varying types of error indicates that this method of determining individual fractionation factors for N^α and N^β is robust even at high levels of error.

3.3 Application of the Expanded Rayleigh model to previously published isotopic data

As a proof-of-concept example, we applied the standard Rayleigh model and the Expanded Rayleigh model to previously published experimental data and compared the results. We chose data collected for N₂O biosynthesis by a culture of M. trichosporium (Methylocystis sp.) grown aerobically in the presence of NH₂OH (Sutka et al., 2006). Methylocystis species encode the gene for hydroxylamine oxidoreductase (HAO), an enzyme that not only converts NH₂OH to NO during nitrification, but also produces N₂O under certain conditions (Yamazaki et al., 2014). As Methylocystis species lack the gene for cytochrome P460 (which has also been shown to convert NH₂OH to N₂O) (Caranto et al., 2016), the isotopic fractionation observed in this experiment is assumed to be due to only one process: HAO-catalyzed conversion of NH₂OH to N₂O. Indeed, the ${\mathit{\epsilon }}_{\text{N-bulk}}$ value for N₂O produced in this experiment (5.3 $\pm$ 0.4 ‰ for replicate B, KIE ¹⁵N${}^{\mathrm{bulk}}=$ 0.9947 $\pm$ 0.0004; see Table S15; Sutka et al., 2006) is similar to the ${\mathit{\epsilon }}_{\text{N-bulk}}$ value of 2.0 ‰ (KIE ¹⁵N${}^{\mathrm{bulk}}=$ 0.9980) reported for purified HAO from Nitrosomonas europaea (Yamazaki et al., 2014). The differences between ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β (SP) are also very similar for the two experiments, indicating that HAO was the primary source of N₂O production in the cultures of M. trichosporium, meaning the Rayleigh model could be applied to this system.

Application of the standard Rayleigh model (Eq. 2) to ${\mathit{\delta }}^{\mathrm{15}}$N^bulk via linear regression of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk against [$-f\ln f/(\mathrm{1}-f)$] indicated that ${\mathit{\epsilon }}_{\text{N-bulk}}$ was 5.3 $\pm$ 0.4 ‰ and KIE ¹⁵N^bulk was 0.9947 $\pm$ 0.0004 (inverse isotope effect) (Table 5) for N₂O production by M. trichosporium. Similarly, application of the standard Rayleigh model to ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β produced a KIE ¹⁵N^α value of 0.9952 $\pm$ 0.0015 and a KIE ¹⁵N^β value of 0.9942 $\pm$ 0.0021. Thus, the standard Rayleigh model predicts an inverse isotope effect for both N^α and N^β, which theoretically means that both positions are enriched in ¹⁵N relative to the substrate. However, all of the ${\mathit{\delta }}^{\mathrm{15}}$N^β values ($-$11.7 ‰ to $-$13.3 ‰) are more negative than the measured ${\mathit{\delta }}^{\mathrm{15}}$N^s0 value of $-$2.3 ‰ (Fig. 4 and Table S14), indicating that ¹⁵N is depleted at the $\mathit{\beta }$ position and that ¹⁴N is preferred (normal isotope effect). The ${\mathit{\delta }}^{\mathrm{15}}$N^α values (21.1 ‰ to 22.7 ‰), in contrast, are more positive than ${\mathit{\delta }}^{\mathrm{15}}$N^s0, indicating that there is an inverse isotope effect at this position. The Expanded Rayleigh model produces position-specific KIEs that match this qualitative assessment: KIE ¹⁵N${}^{\mathit{\alpha }}=$ 0.9779 $\pm$ 0.0004 (inverse), and KIE ¹⁵N${}^{\mathit{\beta }}=$ 1.0121–1.0122 $\pm$ 0.0004 (normal) (Table 5). Additionally, the values of $\mathit{\rho }$ and $\mathit{\tau }$ do not vary significantly with $f$ (Table S15), indicating that fractionation at the $\mathit{\alpha }$ and $\mathit{\beta }$ positions appears to be constant. Comparison of the RMSE values for the standard (linear) Rayleigh model and the nonlinear portion of the Expanded Rayleigh model also indicates that both models have a similar goodness of fit. Thus, the Expanded Rayleigh model can be applied to experimental data to accurately assess position-specific KIEs.

Table 5

Comparison of standard Rayleigh and Expanded Rayleigh KIE values $\pm$ standard error for N₂O produced from NH₂OH by an axenic culture of M. trichosporium (Methylocystis sp.). Values were calculated using isotopic data previously published for M. trichosporium replicate B (Sutka et al., 2006). n/a – not applicable.

Model	Standard Rayleigh	Nonlinear RMSEb	KIE 15Nbulkc	KIE 15Nαc	KIE 15Nβc
	(linear) RMSEa
Standard Rayleighd	0.07	n/a	0.9947 $\pm$ 4 $\times$ 10−4	0.9952 $\pm$ 0.0015	0.9942 $\pm$ 0.0021
Expanded Rayleigh 1e	0.07	0.16	0.9947 $\pm$ 4 $\times$ 10−4	0.9779 $\pm$ 4 $\times$ 10−4	1.0121 $\pm$ 4 $\times$ 10−4
Expanded Rayleigh 2e	0.07	0.16	0.9947 $\pm$ 4 $\times$ 10−4	0.9779 $\pm$ 4 $\times$ 10−4	1.0122 $\pm$ 4 $\times$ 10−4

^a RMSE (root mean square error) for the standard Rayleigh model (Eq. 2) with linear regression of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk against [$-f\ln f/(\mathrm{1}-f)$]. ^b RMSE for nonlinear model 1 or 2 (Eq. 29 or Eq. 30). ^c Average value $\pm$ standard error. ^d KIE values were calculated from ${\mathit{\epsilon }}_{\text{N-bulk}}$, ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}$, or ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$ values obtained via linear regression of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk, ${\mathit{\delta }}^{\mathrm{15}}$N^α, or ${\mathit{\delta }}^{\mathrm{15}}$N^β against [$-f\ln f/(\mathrm{1}-f)$]. ^e For the Expanded Rayleigh model, bulk values (${\mathit{\alpha }}_{\text{N-bulk}}$, ${\mathit{\epsilon }}_{\text{N-bulk}}$, and KIE ¹⁵N^bulk) were determined with the standard Rayleigh approach, $\mathit{\rho }$ was determined via nonlinear regression, and $\mathit{\tau }$ was determined by averaging ¹⁴N${}^{\mathit{\alpha }}{/}^{\mathrm{14}}$N^bulk for every step of the reaction. Then ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$ were calculated with Eq. (21) or Eq. (22) and converted to KIE values using Eq. (10). The only difference between Expanded Rayleigh models 1 (Eq. 29) and 2 (Eq. 30) is which $\mathit{\delta }$ value (${\mathit{\delta }}^{\mathrm{15}}$N^α or ${\mathit{\delta }}^{\mathrm{15}}$N^β) is substituted with a $\mathit{\rho }$-containing expression in the nonlinear model.

We also compared the results of applying the standard and Expanded Rayleigh models to previously published isotopic data on NO reduction to N₂O by purified H. capsulatum P450 NOR (Yang et al., 2014). In this experiment, ${\mathit{\delta }}^{\mathrm{15}}$N^bulk varied linearly with [$-f\ln f/(\mathrm{1}-f)$] ($R^{\mathrm{2}}=$ 0.68–0.94), and averaging the slopes of the standard Rayleigh plots for each individual replicate indicated that ${\mathit{\epsilon }}_{\text{N-bulk}}$ was equal to 9.10 ‰ $\pm$ 1.42 and KIE ¹⁵N^bulk was equal to 0.9910 $\pm$ 0.0014 (Tables 6 and S17). In other words, on average, N₂O became depleted in ¹⁵N (i.e., ${\mathit{\delta }}^{\mathrm{15}}$N^bulk decreased) as the reaction proceeded (Fig. 5 and Table S16). However, the N^α position became more enriched in ¹⁵N as the reaction progressed (i.e., ${\mathit{\delta }}^{\mathrm{15}}$N^α increased as $f$ decreased). Consequently, $\mathit{\rho }$ increased (and $\mathit{\tau }$ decreased) over the course of the reaction (Table S17), and the observed values of ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$ also varied with $f$, in contrast to Assumption 6 (above). Thus, the Expanded Rayleigh model, which is based on the assumption that ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$ remain constant over the course of the reaction, does not fit these data optimally. Therefore, using nonlinear regression (Eq. 29 or Eq. 30) to determine $\mathit{\rho }$ is not appropriate in this case. Using a modified version of the Expanded Rayleigh model, however, we calculated $\mathit{\rho }$ for each individual time point and estimated apparent KIE ¹⁵N^α and KIE ¹⁵N^β values for the beginning and end of the observed range of $f$ values (see the Methods section). Our estimated KIE ¹⁵N^α values ranged from 0.9823 $\pm$ 0.0016 (earlier in the reaction) to 0.9781 $\pm$ 0.0016 (later in the reaction), indicating that the N^α position is subject to an inverse isotope effect throughout the reaction. On the other hand, the KIE ¹⁵N^β estimates ranged from 0.9998 $\pm$ 0.0015 to 1.0041 $\pm$ 0.0013, suggesting that there is essentially no apparent isotope effect at the $\mathit{\beta }$ position, especially early in the reaction. These results are consistent with a qualitative assessment of the $\mathit{\delta }$ values. For example, the ${\mathit{\delta }}^{\mathrm{15}}$N^α values ($-$25.8 ‰ to $-$31.0 ‰) are more positive than the calculated ${\mathit{\delta }}^{\mathrm{15}}$N^s0values (i.e., the intercept of the standard Rayleigh plots for ${\mathit{\delta }}^{\mathrm{15}}$N^bulk), which range from $-$44.7 ‰ to $-$48.3 ‰. This indicates that the $\mathit{\alpha }$ position is enriched in ¹⁵N relative to the substrate. On the other hand, the ${\mathit{\delta }}^{\mathrm{15}}$N^β values near the start of the reaction ($-$45.2 ‰ to $-$48.3 ‰) are similar to the calculated ${\mathit{\delta }}^{\mathrm{15}}$N^s0 values, suggesting that there is very little isotopic fractionation at the $\mathit{\beta }$ position.

In contrast to the Expanded Rayleigh model estimates, the standard Rayleigh calculations yield a normal isotope effect for N^α (KIE ¹⁵N${}^{\mathit{\alpha }}=$ 1.0127 $\pm$ 0.0030) and an inverse isotope effect for N^β (KIE ¹⁵N${}^{\mathit{\beta }}=$ 0.9694 $\pm$ 0.0022). (Note that applying the standard Rayleigh model to individual observations and averaging KIEs for the early and later parts of the reaction yields similar results to applying the standard Rayleigh model to the entire reaction (Tables 6 and S17; see the Supplement for details). While it is possible that N^β exhibits an apparent inverse isotope effect near the beginning of the reaction, the magnitude of this effect is likely much smaller (i.e., early KIE ¹⁵N^β is much closer to 1) than the standard Rayleigh calculations suggest. More strikingly, the standard Rayleigh prediction that N^α is subject to a normal isotope effect (¹⁴N preferred) is clearly incorrect, as the $\mathit{\alpha }$ position is enriched in ¹⁵N relative to the initial substrate. Moreover, as ¹⁵N is depleted from the substrate pool, the fraction of ¹⁵N incorporated into the $\mathit{\alpha }$ position increases, indicating that the preference for ¹⁵N at N^α increases (becomes more inverse) as the reaction proceeds. Overall, analysis of the P450 NOR data demonstrates that even for reactions where the apparent position-specific KIEs vary, the Expanded Rayleigh model produces considerably more accurate KIE estimates than the standard Rayleigh model.

4 Discussion

4.1 Proper application of the Expanded Rayleigh model

The Rayleigh model has been used to study kinetic isotope effects for decades and has been adapted many different times for use in specific scenarios. However, to our knowledge, this is the first examination of how the Rayleigh model applies to a process where two heavy atoms from the same substrate pool (e.g., NO or NH₂OH) are incorporated into non-equivalent positions in the product (e.g., N₂O). Our simulations of different N₂O-producing scenarios demonstrate that the original Rayleigh model should not be applied to the individual atoms in such cases as this produces inaccurate results. However, the new Expanded Rayleigh model presented here can be used to accurately determine individual kinetic isotope effects. This expanded model allows for a more granular analysis than would be possible from determining only bulk kinetic isotope effects and SP.

In practice, when can the Expanded Rayleigh model be used? To develop the Expanded Rayleigh model, we employed the same assumptions that the original model is founded on. Therefore, any situations where the original Rayleigh model applies will also likely be suitable for the Expanded Rayleigh model. More specifically, the Expanded Rayleigh model will likely produce accurate results for any closed system where N₂O is produced from a substrate pool with ¹⁵N at or near natural abundance, as long as the rate-determining step for each atom remains constant. Thus, the Expanded Rayleigh model is ideally suited for studies of isolated N₂O-synthesizing enzymes, as this experimental setup provides a closed system where N₂O is produced by a single catalytic process.

One method for confirming that the rate-determining steps of N₂O synthesis are constant is to determine if ${\mathit{\epsilon }}_{\text{N-bulk}}$ and SP are constant for a particular dataset. If enzymatic catalysis is rate-limiting throughout a time course experiment, ${\mathit{\epsilon }}_{\text{N-bulk}}$ will be constant, and the slope of a standard Rayleigh plot should be linear. In contrast, if substrate diffusion (or some other process) becomes rate-limiting, curvilinear behavior may be exhibited in a standard Rayleigh plot, as has been observed for multi-step N₂O production in microbial cultures (Sutka et al., 2008; Haslun et al., 2018). Of course, to use the Expanded Rayleigh model, ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$ must be constant as well. To verify that the individual fractionation factors are constant, one could plot ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β against [$-f\ln f/(\mathrm{1}-f)$] (i.e., create standard Rayleigh plots for ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β) and confirm that these plots are linear and roughly parallel. Alternatively, SP (the difference of ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β; Eq. 32) can be monitored over time. The lack of a significant trend in SP values is a reasonably good indication that ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$ are constant. In studies of microbial cultures where N₂O is predominantly produced by a single type of enzyme (e.g., bacterial NOR, fungal NOR, or HAO), SP generally remains constant as substrate is consumed (Toyoda et al., 2005; Sutka et al., 2006; Sutka et al., 2008; Haslun et al., 2018). The fact that SP is constant for a variety of N₂O production pathways in a microbial cell suggests that enzyme-specific fractionation at N^α and N^β typically remains constant, meaning the Expanded Rayleigh model can likely be accurately applied to a variety of N₂O-producing enzymes.

Table 6

Comparison of standard Rayleigh and Expanded Rayleigh KIE values $\pm$ standard error for N₂O production from NO by purified Histoplasma capsulatum (fungal) P450 NOR (calculated using previously published isotopic data; Yang et al., 2014).

Model	Extent of reaction		KIE 15Nbulka	KIE 15Nαa	KIE 15Nβa
	(range of $f$)
Standard Rayleighb	All	0.42–0.87	0.9910 $\pm$ 0.0014	1.0127 $\pm$ 0.0030	0.9694 $\pm$ 0.0022
Standard Rayleighc	Early	0.77–0.81	0.9908 $\pm$ 0.0013	1.0130 $\pm$ 0.0031	0.9687 $\pm$ 0.0030
Standard Rayleighc	Late	0.47–0.52	0.9909 $\pm$ 0.0015	1.0121 $\pm$ 0.0029	0.9698 $\pm$ 0.0018
Expanded Rayleighd	Early	0.77–0.81	0.9910 $\pm$ 0.0014	0.9823 $\pm$ 0.0016	0.9998 $\pm$ 0.0015
Expanded Rayleighd	Late	0.47–0.52	0.9910 $\pm$ 0.0014	0.9781 $\pm$ 0.0016	1.0041 $\pm$ 0.0013

^a Average value $\pm$ standard deviation. ^b KIE values were calculated from ${\mathit{\epsilon }}_{\text{N-bulk}}$, ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}$, or ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$ values obtained via linear regression of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk, ${\mathit{\delta }}^{\mathrm{15}}$N^α, or ${\mathit{\delta }}^{\mathrm{15}}$N^β against [$-f\ln f/(\mathrm{1}-f)$]. The standard Rayleigh model values presented here differ slightly from the previously published values (Yang et al., 2014) due to our exclusion of the earliest observation(s) from each replicate (i.e., observations with the highest values of $f$ were excluded). ^c For the standard Rayleigh model applied to individual observations, ${\mathit{\epsilon }}_{\text{N-bulk}}$, ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}$, or ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$ values were determined using Eq. (S24); the $y$ intercept listed in that equation corresponds to the $y$ intercept of ${\mathit{\delta }}^{\mathrm{15}}$N^bulk, ${\mathit{\delta }}^{\mathrm{15}}$N^α, or ${\mathit{\delta }}^{\mathrm{15}}$N^β against [$-f\ln f/(\mathrm{1}-f)$] (determined by linear regression of the data from each replicate). KIE values for six (early) or seven (late) individual observations were pooled and averaged. ^d For the Expanded Rayleigh model applied to individual observations, bulk values (${\mathit{\alpha }}_{\text{N-bulk}}$, ${\mathit{\epsilon }}_{\text{N-bulk}}$, and KIE ¹⁵N^bulk) were determined with the standard Rayleigh approach. $\mathit{\rho }$ was calculated for each observation using Eq. (23) ($\mathit{\rho }=$ ¹⁵N${}^{\mathit{\alpha }}{/}^{\mathrm{15}}$N^bulk), and $\mathit{\tau }$ was determined for every step of the reaction using Eq. (24) ($\mathit{\tau }{=}^{\mathrm{14}}$N${}^{\mathit{\alpha }}{/}^{\mathrm{14}}$N^bulk). Then ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\alpha }}_{\mathrm{N}-\mathit{\beta }}$ were calculated for each individual observation with Eq. (21) or Eq. (22) and converted to KIE values using Eq. (10). KIE ¹⁵N^α and KIE ¹⁵N^β values for six (early) or seven (late) individual observations were pooled and averaged.

At present, there is a scarcity of studies reporting N₂O fractionation by purified enzymes for multiple values of $f$, limiting our ability to test the Expanded Rayleigh model with published experimentally derived datasets. We were, however, able to apply the Expanded Rayleigh model to previously published data on M. trichosporium cultures (Sutka et al., 2006), whose HAO-catalyzed oxidation of NH₂OH was the primary source of N₂O and displayed constant values of ${\mathit{\epsilon }}_{\text{N-bulk}}$ and SP. Thus, we were able to calculate observed values of KIE ¹⁵N^α and KIE ¹⁵N^β that serve as reasonable estimates of HAO-specific isotope effects (Table 5). Additionally, we have used the Expanded Rayleigh model to re-evaluate isotopic data our group has collected for purified fungal P450 NOR (Table 6) (Yang et al., 2014). The value of ${\mathit{\epsilon }}_{\text{N-bulk}}$ for P450 NOR appears to be constant, but SP increased by 14 ‰ as the fraction of NO reduced changed from 10 % to 50 % (Yang et al., 2014). Given the proposed P450 NOR reaction mechanism, the observed isotope effect for the first NO molecule is influenced by both an equilibrium isotope effect related to NO binding and an intrinsic KIE associated with catalysis (Yang et al., 2014). Conversely, the rate-determining steps for the second NO molecule are likely only associated with catalysis. Thus, the isotopic fractionation for the first nitrogen atom to be incorporated will be more sensitive to changes in NO concentration as the reaction progresses than the second N atom, leading to SP values that are not constant. Even in these circumstances, however, the Expanded Rayleigh model estimates for KIE ¹⁵N^α and KIE ¹⁵N^β are a significant improvement on the previously published standard Rayleigh estimates, which incorrectly predict a normal isotope effect at the $\mathit{\alpha }$ position. Overall, the Expanded Rayleigh model fits the available data on isotopic fractionation during N₂O biosynthesis reasonably well, although we recommend collecting isotopic data across a range of $f$ values to confirm that SP (and therefore KIE ¹⁵N^α and KIE ¹⁵N^β) remains constant when using this model.

As outlined above (in the Assumptions section), the Expanded Rayleigh model is most accurate for reactions where the isotopic fractionation factors (${\mathit{\epsilon }}_{\text{N-bulk}}$, ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}$, and ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$) and ${\mathit{\delta }}^{\mathrm{15}}$N^s0 values are not too extreme. This stipulation, which also applies to the standard Rayleigh equation, is necessary in part because an approximation introduced by Mariotti and colleagues (Mariotti et al., 1981) decreases in accuracy the farther ${\mathit{\delta }}^{\mathrm{15}}$N^s0 and ${\mathit{\delta }}^{\mathrm{15}}$N^s are from 0 ‰. Additionally, extremely high positive values of ${\mathit{\delta }}^{\mathrm{15}}$N^s0 or ${\mathit{\epsilon }}_{\text{N-bulk}}$ would result in large positive product $\mathit{\delta }$ values, corresponding to an increased abundance of ¹⁵N in the product. As both the nonlinear Rayleigh model and the underlying calculation of ${\mathit{\delta }}^{\mathrm{15}}$N^β values depend on the assumption that the amount of ¹⁵N in N₂O is close to the natural abundance, the accuracy of the Expanded Rayleigh model would decrease in such circumstances. Thus, applying the Expanded Rayleigh model to experiments with spiked ¹⁵N is not advised. Simulations of $\mathit{\delta }$ values calculated without using Mariotti's approximation indicate that for absolute values of ${\mathit{\epsilon }}_{\text{N-bulk}}$ or ${\mathit{\delta }}^{\mathrm{15}}$N^s0 less than or equal to 50 ‰, the level of error introduced by Mariotti's approximation is similar to the expected level of experimental error (i.e., absolute relative difference in KIE values of $\sim$ 0.001) (Tables S1–2).

The magnitude of ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\alpha }}$ and ${\mathit{\epsilon }}_{\mathrm{N}-\mathit{\beta }}$ must also be considered when applying the Expanded Rayleigh model. In theory, if isotopic enrichment at N^α or N^β is very large, the assumption that essentially equal amounts of ¹⁴N are apportioned to each position (i.e., that $\mathit{\tau }\approx \mathrm{0.5}$) would be violated, and Eq. (15) would no longer hold. However, violating this assumption would require unusually large heavy-atom kinetic isotope effects or a spiked sample. As outlined in Table S3, for a simulated reaction where ${\mathit{\epsilon }}_{\text{N-bulk}}$ is set at $-$20 ‰ (KIE ¹⁵N${}^{\mathrm{bulk}}=$ 1.0204) and ${\mathit{\delta }}^{\mathrm{15}}$N^s0 is 0 ‰, individual KIE values would have to be significantly lower than 0.96 or greater than 1.08 to introduce error greater than analytical error. These KIE ¹⁵N^α or KIE ¹⁵N^β values correspond to $\mathit{\rho }$ values less than 0.47 or greater than 0.53. In practice, our group has found that for two different types of N₂O-producing enzymes, $\mathit{\rho }$ is nowhere near these extremes. For example, for purified P450 NOR (Yang et al., 2014) and for HAO-catalyzed NH₂OH oxidation in an axenic culture (Sutka et al., 2006), $\mathit{\rho }$ does not exceed 0.509 (Tables S15 and S17). Overall, we anticipate that the Expanded Rayleigh model will be applicable to a wide variety of N₂O biosynthesis reactions.

4.2 Other practical considerations

The Expanded Rayleigh model presented here does not require measurement of ${\mathit{\delta }}^{\mathrm{15}}$N^s0. This is advantageous for N₂O synthesis reactions where NO is the substrate because directly measuring an accurate ${\mathit{\delta }}^{\mathrm{15}}$N value for NO is quite challenging due to the reactivity of NO, the lack of isotopic standards for NO, and the fact that ¹⁴N¹⁶O has the same molecular weight as ¹⁵N¹⁵N. However, ${\mathit{\delta }}^{\mathrm{15}}$N^s0 can be directly measured when NH₂OH is used as the substrate for N₂O production. Determining ${\mathit{\delta }}^{\mathrm{15}}$N^s0 experimentally may be useful in validating application of the standard Rayleigh model, as the intercept of the standard Rayleigh model is expected to be equal to ${\mathit{\delta }}^{\mathrm{15}}$N^s0. For example, we determined that the intercept of the standard Rayleigh plot for N₂O produced by NH₂OH oxidation by M. trichosporium is 0.5 $\pm$ 0.3 ‰, which is similar to the measured ${\mathit{\delta }}^{\mathrm{15}}$N^s0 value of $-$2.3 ‰ (Sutka et al., 2006), confirming that the standard Rayleigh model is applicable for this dataset. Thus, when feasible, measuring ${\mathit{\delta }}^{\mathrm{15}}$N^s0 would complement the expanded Rayleigh approach.

The nonlinear portion of the Expanded Rayleigh model (used to determine $\mathit{\rho }$) can be written in two forms, nonlinear model 1 (Eq. 29) and nonlinear model 2 (Eq. 30). Both models are based on the assumption that ${\mathit{\delta }}^{\mathrm{15}}$N^bulk is equal to the average of ${\mathit{\delta }}^{\mathrm{15}}$N^α and ${\mathit{\delta }}^{\mathrm{15}}$N^β (Eq. 15). Therefore, both models are expected to produce very similar results as long as that assumption holds. Indeed, as we have demonstrated here, both models yield nearly identical values of $\mathit{\rho }$, KIE ¹⁵N^α, and KIE ¹⁵N^β for simulated datasets (Figs. 6 and S5–8; Tables S9–13), and we have obtained similar results for the available experimental data on N₂O production from NH₂OH (Table 5). Thus, either version of the Expanded Rayleigh model may be used. If there is any doubt about which version to use, we recommend selecting the model that has a lower RMSE value (and thus fits the data better).

5 Conclusions

In this paper, we presented the Expanded Rayleigh model, a novel adaptation of the Rayleigh distillation equation that describes the position-specific isotopic enrichment that occurs at N^α and N^β during N₂O synthesis. Using simulated datasets representing multiple different combinations of KIEs, we demonstrated that the Expanded Rayleigh model accurately recapitulates KIE ¹⁵N^α and KIE ¹⁵N^β for a wide variety of scenarios. Our simulations also demonstrate that this new model is robust even when applied to skewed data and/or data with a high level of error. Additionally, we have shown the Expanded Rayleigh model fits the experimentally measured isotopic data for N₂O production from NH₂OH quite well (Sutka et al., 2006) and provides significantly improved estimates for the position-specific KIEs for P450 NOR-catalyzed N₂O synthesis. Thus, the Expanded Rayleigh model provides a reliable method for quantifying isotopic fractionation at N^α and N^β, and it promises to be a valuable tool for experimentally probing the catalytic mechanisms of N₂O-synthesizing enzymes. Finally, we note that the Expanded Rayleigh model will also likely be appropriate for the analysis of other reactions where two atoms from the same substrate pool are incorporated into distinct, non-exchangeable positions in the product.

Code and data availability

All code used in this manuscript is available at 10.5281/zenodo.13931602 (Rivett, 2024a). The results for each of the 1000 simulations with error (derived from Datasets 1–5) can be found at 10.5281/zenodo.13323399 (Rivett, 2024b).

The supplement related to this article is available online at: https://doi.org/10.5194/bg-21-4549-2024-supplement.

Author contributions

EDR, NEO, and ELH conceived this study. EDR developed the model and analyzed the data with contributions from WM, NEO, and ELH. WM provided critical input for statistical analysis and coding. EDR prepared the original draft of the manuscript and figures, and all authors participated in editing and revision of the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

We thank Julius Campeciño and Helen Kreuzer for a careful reading of a draft of this manuscript and for their insightful comments and questions. We also thank Julius Campeciño for helpful discussions during initial development of the Expanded Rayleigh model. We also acknowledge an anonymous reviewer whose comments first alerted us to the complications associated with position-specific application of the Rayleigh model.

Financial support

This material is based upon work supported by the Great Lakes Bioenergy Research Center, US Department of Energy, Office of Science, Biological and Environmental Research Program under award number DE-SC0018409.

Review statement

This paper was edited by Perran Cook and reviewed by one anonymous referee.