Metabolic control analysis of biogeochemical systems

Many reactive systems comprise several processes operating at different spatial and temporal scales, such as cellular metabolism and regulation, enzyme kinetics, microbial population dynamics, abiotic chemical transformations, hydrodynamic and diffusive transport, downward migration of matter in water and sediment columns, tectonic uplift and erosion. Examples of such systems include local or global nitrogen cycling¹, the ocean carbon pump², microbial mats³, material and energy fluxes in food webs⁴, bioreactors, and wastewater treatment plants⁵. Understanding the sensitivity of a system’s emergent properties to each of its constituent processes, and identifying the subset of processes that are truly relevant for predictions, is important for multiple reasons. First, parameterizing each of the aforementioned processes in models requires considerable effort; for example, determining microbial metabolic and growth kinetics typically involves laborious incubation experiments^6,7. In some microbial systems, physical mixing across space can be the main rate-limiting process, effectively eliminating the dependency on microbial kinetic and population-dynamic parameters^8,9. Second, models that account for processes operating at vastly different spatial or temporal scales tend to be computationally expensive^10,11,12, and thus shortcutting some secondary processes could allow for considerable optimization. Third, efficiently steering a system in a specific desired direction, for example, to boost marine productivity or increase wastewater treatment rates, necessitates a determination of the critical “turning knobs”¹³. Flexible and practical mathematical frameworks for performing such analyses are thus highly desirable.

For simple cellular biochemical reaction networks at steady state or quasi-steady-state, metabolic control analysis (MCA) has become a standard mathematical framework for examining the controls that individual enzyme activities exert on fluxes and substrate concentrations^14,15,16. Partial derivatives of the network’s responses, including substrate concentrations and fluxes, with respect to enzyme activities are used to define local sensitivity coefficients called concentration control coefficients (CCCs) and flux control coefficients (FCCs), respectively. Due to the particular structure of reaction networks, MCA yields mathematical theorems facilitating the interpretation of FCCs and CCCs, and provides practical guidelines for experimentally determining these coefficients. A notable conceptual impact of MCA has been the realization that many biochemical reaction networks are not governed by a single rate-limiting step, and that their control is instead distributed over multiple steps¹⁷. MCA has been instrumental in studies of cell metabolism, metabolic engineering, and drug development, to name a few applications^18,19. To date, MCA has not been generalized to reactive systems with explicit spatial structure or involving processes other than chemical reactions, such as reaction-advection-diffusion models for carbon cycling in sediment columns, microbial communities distributed across a water column, global ocean circulation-biogeochemical models and bioreactors with engineered control loops. Such a generalization would yield novel insights of similar impact as those achieved by MCA in biochemistry. For example, a generalized MCA would facilitate the examination of controls exerted by processes other than enzymatic reactions, such as molecular diffusion across space, microbial growth, and death, hydrodynamic mixing and sedimentation in the ocean, tectonic uplift in Earth-system models, or hydraulic renewal in wastewater treatment plants.

Here I show how MCA can be elegantly generalized to a much wider class of reactive systems than traditionally considered, including all of the examples mentioned above. In particular, I identify a simple set of conditions, which are widely encountered in reactive systems and which are sufficient for the applicability of MCA. In a nutshell, the most crucial conditions are the existence of a subset of focal parameters, whose uniform rescaling leaves the system’s steady state unchanged but rescales fluxes by the same factor. As I will show, these parameters are analogous to enzyme activities in biochemical networks analyzed through classical MCA, and a rescaling of these parameters can often be interpreted as a simple rescaling of time. I further define the FCCs and CCCs in a more general form than commonly done for biochemical networks, prove the validity of two important theorems known as the flux and concentration control summation theorems, and provide various examples of models amenable to MCA. I then demonstrate this generalized MCA in detail for two marine systems, each described by a reaction-advection-diffusion partial differential equation model: The sulfate methane transition zone in Black Sea sediments^20,21, where I examine the controls of various processes on anaerobic methane oxidation rates, and a seasonal oxygen minimum zone at quasi-steady state in Saanich Inlet^22,23,24, where I examine the controls on microbially driven fixed nitrogen loss.

Generalizing MCA

I henceforth provide a formal specification of the systems to which MCA shall be generalized, definitions of FCCs and CCCs for such systems, and simple proofs of the famous flux and concentration control summation theorems. As a recurring illustrative example, I will consider a general reaction-advection-diffusion partial differential equation (PDE) model on a finite 1-dimensional interval:

where z represents the single spatial coordinate within the interval ({{{mathcal{Z}}}}=[a,b]), the dynamic variables U₁, . . , U_m are scalar fields on ({{{mathcal{Z}}}}) defining the system’s state at any moment in time, D₁, . . , D_m are given location-dependent diffusivities, v₁, . . , v_m are given location-dependent advection speeds and S₁, . . , S_m are source terms (e.g., accounting for chemical or biological processes). For example, z may represent depth along a water or sediment column, and each U_j may represent the concentration of a specific substance or the population density of a microbial species. While the above model is a priori dynamical, I will henceforth focus on its steady states, i.e., the case where ∂U_j/∂t = 0. For illustration, at the first boundary I consider fixed-value conditions:

whereas at the second boundary I consider fixed-flux conditions:

where A₁, . . , A_m and B₁, . . , B_m are constants. Other types of boundary conditions, including more general Robin boundary conditions²⁵, are also possible. Additional examples are given throughout the article and in Supplementary Discussion 1. In particular, Supplementary Discussion 1.1 elaborates on how the general concepts presented here relate to classical MCA, while Supplementary Discussion 1.3 exemplifies these concepts on a 3-dimensional carbon cycle PDE model overlaid on an ocean general circulation model.

State space: I consider dynamical systems at (quasi-)steady state, whose state at any point in time is described by a vector-valued function U defined on some arbitrary domain ({{{mathcal{Z}}}}), i.e., ({{{bf{U}}}}:{{{mathcal{Z}}}}to {{mathbb{R}}}^{m}). Let ({{{mathcal{U}}}}) denote the set of all possible system states. Here, m represents the number of dynamic physical/biological variables, such as chemical concentrations or biological species densities. For a well-mixed system described by an ordinary differential equation model, such as biochemical reaction networks considered in classical MCA¹⁵, ({{{mathcal{Z}}}}) is a single point and thus U is just a vector of m numbers, whereas for an ordinary differential equation box-model with multiple connected boxes, ({{{mathcal{Z}}}}) is an enumeration of the boxes and thus U comprises m numbers for each of the boxes. In the case of a diffusion-advection PDE model in 3 spatial dimensions, ({{{mathcal{Z}}}}) represents the 3 spatial axes (i.e., ({{{mathcal{Z}}}}) is a subset of ({{mathbb{R}}}^{3})), U represents m scalar fields defined at each point in space, and ({{{mathcal{U}}}}) represents the space of acceptable fields (for example, the set of twice-continuously differentiable functions ({{{mathcal{Z}}}}to {{mathbb{R}}}^{m})²⁶). For the PDE example given in Eq. (1), ({{{mathcal{Z}}}}) is simply the interval [a, b] and each U represents m curves on [a, b].

Henceforth, let ({{{bf{p}}}}in {{mathbb{R}}}^{n}) denote n focal model parameters of interest. Note that p does not need to comprise all conceivable model parameters, but rather represents a focal subset of n parameters with respect to which we will be examining the system’s sensitivity. All other model parameters not included in p are assumed fixed throughout. In the PDE model in Eq. (1), p could for example include the assumed advection velocities, the imposed boundary fluxes B_j, and so on. Later, I will discuss specific conditions that p must satisfy in order for MCA to be applicable.

Definition of steady state: Assume that for any given value of p, the steady state is uniquely defined as the solution to the following set of abstract equations:

where the index ω takes values in some arbitrary set Ω, and each F_ω is a scalar-valued function that imposes a condition on U, i.e., ({F}_{omega }:{{{mathcal{U}}}}times {{mathbb{R}}}^{n}to {mathbb{R}}). The form of each F_ω can be quite broad, including algebraic, differential and even integral equations. Similarly, Ω may in principle be arbitrary, although in practice it will be closely related to the domain ({{{mathcal{Z}}}}) and the number of variables m. For example, in the case of a single-box ordinary differential equation model, Ω would be {1, . . , m} and F₁, . . , F_m would constitute m algebraic equations that fully determine the steady state. In the case of a PDE model on a finite 3-dimensional domain, Ω would be ({{{mathcal{Z}}}}times {1,..,m}), with F_ω encoding the PDE in the interior of ({{{mathcal{Z}}}}) as well as the boundary conditions at the boundary of ({{{mathcal{Z}}}}). For the PDE example in Eq. (1), Ω is [a, b] × {1, . . , m}, with steady-state conditions being:

where for simplicity I omitted writing out the focal parameters p in the various terms above. For any given value of p, let (overline{{{{bf{U}}}}}({{{bf{p}}}})in {{{mathcal{U}}}}) be the corresponding steady state, i.e. satisfying

for all ω ∈ Ω.

Scaling invariance of steady state: We are now ready to specify the first major condition for MCA. Assume that uniformly rescaling the parameters p by an arbitrary positive constant α > 0 does not change the steady state, in other words:

The above invariance condition is the most critical condition that a system and the focal parameters must satisfy, in order for MCA to be applicable. Importantly, whether Eq. (7) is satisfied or not depends on which parameters are included in p. As an example, in biochemical reaction networks considered by classical MCA¹⁵, p typically represents the activities of enzymes, while other parameters such as stoichiometric coefficients are kept fixed. In these models, reaction rates are assumed to scale linearly with enzyme activities when substrate concentrations are kept constant, and steady-state concentrations of metabolites remain unchanged if all enzyme activities are modified by the same factor while keeping all other parameters fixed (details in Supplementary Discussion 1.1). As illustrated in further examples below, basic familiarity with the system and mathematical intuition are often sufficient for recognizing a proper subset of focal parameters, however, in some cases the task may be more complicated. Below, all control coefficients will be defined with respect to these focal parameters p. It is worth mentioning that condition (7) is equivalent to the condition that each F_ω satisfies:

for all α > 0. It is clear that if F_ω is positively homogeneous in p, that is, satisfies:

for some fixed (kappa in {mathbb{R}}) and for all p, all ({{{bf{U}}}}in {{{mathcal{U}}}}) and all α > 0, then F_ω would necessarily also satisfy condition (8). Note that the exponent κ may in principle depend on ω. I stress that homogeneity of F_ω in p is a stronger condition than condition (8) or the equivalent (7), in other words condition (9) always implies (8) and (7), but not the other way around. Nevertheless, as the examples below will show, homogeneity is often satisfied and easy to recognize.

Concentration control coefficients: We are now ready to define the first set of control coefficients, known in classical MCA as concentration control coefficients (CCCs). Consider a scalar quantity of interest C that can be expressed as a function solely of U and not explicitly depending on p, that is, C is fully determined by U and thus (C:{{{mathcal{U}}}}to {mathbb{R}}). For any given p, define (overline{C},({{{bf{p}}}})) as the corresponding steady state value of C, i.e., (overline{C},({{{bf{p}}}}):= C,(overline{{{{bf{U}}}}}({{{bf{p}}}}))). For example, in a biochemical reaction network C may be the concentration of a specific metabolite. In the PDE example in Eq. (1), C may be U_j at a specific location z, i.e., C (U) = U_j(z), or it may be the average value of U_j, i.e., (C,({{{bf{U}}}})=frac{1}{b-a}int_{a}^{b}{U}_{j}(z),dz). To emphasize the analogy to classical MCA, I henceforth refer to C as a concentration regardless of its precise nature. The CCCs of this concentration are defined as:

where i ∈ {1, . . , n} enumerates all focal parameters. In classical MCA, one typically considers the concentrations of all metabolites, and thus one obtains a distinct CCC for each metabolite and for each enzyme. Each μ_i is a unitless number that quantifies the relative change in (overline{C}) under an infinitesimal change in the parameter p_i. In the context of generic sensitivity analysis, such coefficients are also known as normalized local sensitivity coefficients²⁷. Each CCC only encodes a local sensitivity, i.e., the response to small changes, and may (and usually does) depend on p. Recall that we assumed that (overline{{{{bf{U}}}}}) is invariant under a uniform rescaling of p, which means that (overline{C}) is also invariant. By Euler’s homogeneous function theorem²⁸, it follows that:

Dividing Eq. (11) by (overline{C}) yields the famous concentration control summation theorem:

The theorem essentially asserts that the controls of the various focal parameters on any given “concentration” cancel each other out, with some parameters having negative controls and others having positive controls.

Flux control coefficients: We now consider the control that the focal parameters p exert on fluxes. As will become more clear in the definition below and in later examples, the meaning of flux is kept quite broad in this article, and may include the rate of a reaction in a biochemical network, or the net consumption rate of a substrate by multiple reactions, or the vertical flux of a substance across a sediment column, or the cell production rate for a specific microbial species. Specifically, suppose that our flux of interest can be expressed as a function (J:{{{mathcal{U}}}}times {{mathbb{R}}}^{n}to {mathbb{R}}) that depends on U and p and that is positively 1-homogeneous in p, that is, it satisfies:

for all α > 0. The above condition is the second-most critical condition that a system must satisfy for MCA of fluxes in this article (the first one being Eq. (7)). Intuitively, the above homogeneity condition states that uniformly rescaling all focal parameters p by a positive factor, while keeping the system state fixed, should modulate the flux by that same factor. For example, in classical MCA, uniformly altering all enzyme activities by the same factor, while keeping metabolite concentrations fixed, increases each reaction rate by that same factor. In the PDE model in Eq. (1), increasing all diffusion coefficients and advection speeds uniformly by the same factor while keeping U unchanged, would increase flux rates across space by that same factor. In some atypical cases, substantial system expertise may be needed to find a good set of focal parameters p in conjunction with proper flux functions so as to satisfy Eq. (13), in addition to satisfying the scaling invariance condition (7) for (overline{{{{bf{U}}}}}). For any given p, let (overline{J}({{{bf{p}}}})) be the flux at steady state:

We are now ready to define a second set of control coefficients, known in classical MCA as flux control coefficients (FCCs):

where i ∈ {1, . . , n} enumerates all focal parameters. Each ν_i is a unitless number that quantifies the relative change in the steady state flux (overline{J}({{{bf{p}}}})) under an infinitesimal change in the parameter p_i. Note that FCCs are specific to the flux considered, i.e., one would obtain a separate set of FCCs if one were to consider a different flux. Similarly to CCCs, a summation theorem can also be derived for FCCs, as follows. For any p and any α > 0 we have:

in other words, (overline{J}) is positively 1-homogeneous. By Euler’s homogeneous function theorem²⁸ it thus follows that:

Dividing both sides of Eq. (17) by (overline{J}) yields the famous flux control summation theorem:

The theorem asserts that the FCCs of a given flux must sum to 1. In many cases (but not always, e.g.²⁹), the FCCs are non-negative numbers, and can thus be interpreted as fractional contributions to the overall control of a flux. For example, in classical MCA the FCCs for any given reaction rate correspond to the controls exerted by the reaction network’s individual enzymes, which in many cases (e.g., for linear pathways) are all non-negative. In such cases no more than one of the parameters may exert full control with an FCC near 1. At least in classical MCA, typically all FCCs are much lower than 1, i.e., none of the enzymes is of major importance in and of itself, and major changes in flux are only achievable by simultaneously changing multiple enzyme activities^17,29. Whether this situation also turns out to be common in other systems, accessible to the generalized MCA described here, is an open question (e.g., see simulation examples below).

Intuitive interpretation of focal parameters: It should be clear now that the applicability of MCA, including the CCC and FCC summation theorems, requires the recognition of a proper set of parameters, here referred to as focal. In principle, there may exist multiple alternative sets of focal parameters in a model, and this is more likely to occur in more redundantly parameterized models. Typically, one set of focal parameters will stand out as more intuitive, interesting and/or less trivial. In many biogeochemical models, including the two examples discussed in detail below as well as the examples in Supplementary Discussion 1, an intuitive and non-trivial set of focal parameters can be obtained by considering all rate-like parameters such as microbial per-capita death rates, reaction rate coefficients, outgassing rate coefficients, advection velocities, diffusivities, sedimentation rates and so on, as well as the inverses of all time-like parameters such as half-lifetimes of decay processes. Such parameters can usually (but not always) be readily recognized based on their physical units, for example because they explicitly include time or inverse time (e.g., seconds or per-second) or because they can be redefined in such a manner (e.g., Watt = Joules per second). Generally, such parameters determine the time scales of various processes in the system, and ultimately set the “pace” of the system’s dynamics. Multiplying all such parameters by the same factor α thus corresponds to a simple linear rescaling of time, which in and of itself leaves the steady state unchanged (scaling invariance condition in Eq. (7)) but modulates process rates and fluxes by that factor (homogeneity condition in Eq. (13)). Similar time-scaling arguments have been made before for classical MCA^30,31,32. FCCs with respect to such parameters yield insight into the rate-limiting effects that individual processes potentially operating at different time scales (e.g, chemical reactions, physical transport, biological growth) have on energy and material fluxes in a system. Note that while the above serves as a rough intuitive guideline for finding a set of focal parameters, the conditions discussed above (Eqs. (7) and (13)) should always be checked for any candidate set of focal parameters.

Grouping parameters: In some cases, the number of parameters that must be included in p in order to achieve scaling invariance of the steady state (Eq. (7)) and homogeneity of fluxes (Eq. (13)) can be quite large, and the control coefficients associated with each such parameter will tend to be close to zero. For example, in reaction-advection-diffusion models such as the PDE in Eq. (1) diffusivity may vary across space, in which case diffusivity alone comprises multiple parameters — in the extreme case, a separate diffusivity per grid point^33,34. In such cases it can be helpful to combine multiple similar parameters into a single group to evaluate their overall control on the system. In the following I explain how such a grouping might be achieved, and how the control coefficients of a parameter group relate to the control coefficients of the original grouped parameters.

Let p be a set of focal parameters satisfying the MCA assumptions discussed above (Eqs. (7) and (13)), and suppose that we wish to combine the first k parameters, i.e., p₁, . . , p_k, into a single new parameter ({tilde{p}}_{1}). This representative parameter may be set, for example, to the average of the original p₁, . . , p_k, although many other representations are conceivable. We henceforth define ({q}_{i}:= {p}_{i}/{tilde{p}}_{1}) for i = 1, . . , k, so that we can express the original parameters as ({p}_{i}={tilde{p}}_{1}cdot {q}_{i}). We thus obtain a reduced set of n − k + 1 focal parameters, ({tilde{p}}_{1},{tilde{p}}_{2}:= {p}_{k+1},..,{tilde{p}}_{n-k+1}:= {p}_{n}), collectively denoted as a vector (tilde{{{{bf{p}}}}}). Importantly, while the q₁, . . , q_k remain essential for the full model specification, they are henceforth kept fixed and are no longer considered for MCA. Thus, we can ask what happens to the steady state when ({tilde{p}}_{1}) is varied while the q₁, . . , q_k remain fixed. The new focal parameters (tilde{{{{bf{p}}}}}) again satisfy the MCA assumptions (Eqs. (7) and (13)). Indeed, uniformly rescaling (tilde{{{{bf{p}}}}}) by a factor α is equivalent to rescaling the original p by that factor, which by assumption keeps (overline{{{{bf{U}}}}}) unchanged and changes J by the factor α. The CCC of any concentration with respect to ({tilde{p}}_{1}) is:

Thus, a CCC with respect to ({tilde{p}}_{1}) is simply the sum of CCCs with respect to the original represented parameters p₁, . . , p_k. Similarly, the FCC with respect to ({tilde{p}}_{1}) is:

Thus, an FCC with respect to ({tilde{p}}_{1}) is the sum of FCCs with respect to the original represented parameters. Intuitively speaking, the control exerted by ({tilde{p}}_{1}) on a steady state concentration or a steady state flux is the sum of controls that would be exerted individually by the p₁, . . , p_k.

Example: Sulfate-methane transition zone in Black Sea sediments

As a first detailed example, I examine a simple PDE model for a sulfate-methane transition zone (SMTZ) in Black Sea sediments^20,21. SMTZs are widespread zones in marine sediments where methane diffusing upward from deep sediments is oxidized by microbial consortia using downward moving sulfate as a terminal electron acceptor³⁵. The sulfate itself is reduced to hydrogen sulfide, a toxin to many organisms, which accumulates locally but is ultimately consumed by processes on both ends of the SMTZ. SMTZs effectively act as a barrier to methane that might otherwise be released into the water and ultimately the atmosphere. The efficacy of this barrier in principle depends on the kinetics of local chemical reactions involved in the sulfur and methane cycles, as well as rates of diffusive/advective transport across the sediment column. The goal of this section is to examine the controls that these factors exert on methane oxidation rates and sulfide accumulation.

The model is based on insights by Jørgensen et al.²⁰ and Egger et al.²¹, and captures the major microbial and physical processes governing methane and sulfur fluxes in the top 600 cm of sediments. Specifically, the model is a system of 1-dimensional reaction-advection-diffusion PDEs across the sediment column, for the concentrations of m = 6 dynamic metabolites: Methane (CH₄), sulfate (({{{{rm{SO}}}}}_{4}^{2-})), hydrogen sulfide (H₂S), divalent iron (Fe²⁺), iron sulfide (FeS) and pyrite (FeS₂). All metabolites are subject to downward advection due to sedimentation, while dissolved metabolites are additionally transported via molecular diffusion. Molecular diffusivities are adjusted for sediment porosity²¹ and are thus depth-dependent. Metabolites are coupled through 3 chemical reactions: Microbially driven anaerobic oxidation of methane using sulfate (AOM):

formation of iron sulfide (FIS):

and formation of pyrite (FP):

Thus, each metabolite’s concentration U_j (mol ⋅ L⁻¹) satisfies the following PDE:

where z denotes depth, D_j is the depth-dependent diffusivity (cm² ⋅ yr⁻¹) and v_j the depth-dependent downward advection speed (cm ⋅ yr⁻¹) for the jth metabolite, ({mathbb{S}}) is the stoichiometric matrix of the reactions and R_r is the rate of the r-th reaction (mol ⋅ L⁻¹ ⋅ yr⁻¹). While the diffusivities are different for each metabolite, advection speeds only differ between solids and dissolved metabolites; apart from these factors, diffusivities and advection speeds are largely determined by the sedimentation rate, sediment density and the porosity profile (Supplementary Discussion 2). The above model is a priori dynamical, however, we henceforth focus on its steady states, i.e., where ∂U_j/∂t = 0. Reaction kinetics are expressed as mass action laws, i.e., 1st order in each substrate:

where β_r is a constant rate coefficient (henceforth “affinity”) and L_r specifies the substrates for the r-th reaction (e.g., CH₄ and ({{{{rm{SO}}}}}_{4}^{2-}) for methane oxidation). Note that mass action laws can be seen as approximations of common higher-order rate laws when substrate concentrations are low^7,36,37. This model does not track and does not explicitly account for microbial population dynamics. Boundary conditions on both ends, i.e., at a = 0 and b = 600 cm, are in the form of fixed values for all metabolites:

where A_j and B_j are constants. For any given choice of parameters, the corresponding steady state (overline{{{{bf{U}}}}}) was computed by integrating the PDE forward in time until convergence. For additional details, including model parameterization, boundary values and numerical simulation scheme, see Supplementary Discussion 2 and Supplemental Tables S1 and S2. The steady state for default parameters is shown in Fig. 1A–F, with corresponding net production rates of metabolites shown in Fig. 1G–L. Consistent with general expectations^20,21, methane oxidation takes place in a narrow zone at intermediate depth, where methane and sulfate concentrations both converge to zero. Hydrogen sulfide generated by this process either escapes past the sediment-water interface at the top, or is consumed through the formation of pyrite below (Fig. 1E, K). Iron sulfide (Eq. (22)) mostly acts as an intermediate to the eventual formation of pyrite (Eq. (23)).

Steady state conditions and focal parameters: Following the notation introduced in Eq. (5), the system’s steady state conditions are the solutions to the following set of equations:

So far we have not yet specified which of the model parameters should be focal, which is necessary in order to fully apply MCA. As a reminder, focal parameters must at least satisfy the invariance condition for steady states in Eq. (7) or equivalently Eq. (8). If fluxes are to be analyzed, focal parameters must also satisfy the homogeneity condition for fluxes in Eq. (13). Following the earlier guideline on determining a set of focal parameters, one may consider all rate-like parameters in the model, namely the diffusivities D_j, the advection speed v_j and the affinities β_r. For illustration, these are highlighted below using arrows in the first set of steady-state conditions, copied from Eq. (27):

Note that there are no time-like parameters in this model. Uniformly rescaling the D_j, v_j and β_r with a constant factor (while keeping U fixed) rescales the right hand side of Eq. (28) by the same factor, thus satisfying the homogeneity condition in Eq. (9) with κ = 1, for all F_z,j, where a < z < b. Further, rescaling these parameters (again, while keeping U fixed) has no effect on the boundary conditions, i.e., F_a,j(U, αp) = F_a,j(U, p) and F_b,j(U, αp) = F_b,j(U, p), thus satisfying the homogeneity condition in Eq. (9) with κ = 0. Thus, collectively the D_j, v_j and β_r indeed constitute a suitable set of focal parameters. Intuitively, these parameters define the “pace” of the system’s dynamics, and uniformly rescaling these parameters corresponds to a simple rescaling of time, which always preserves steady states. Uniformly rescaling all focal parameters by a factor (while keeping U fixed) also rescales each reaction rate R_r by that same factor, and thus reaction rates satisfy the homogeneity condition for fluxes in Eq. (13). This conclusion remains valid even when considering the mean reaction rates, i.e., averaged over the entire depth range.

I stress that the invariance condition for steady states would no longer hold if one omitted any single one of D_j, v_j or β_r from the set of focal parameters. For example, if one were to only consider the D_j and β_r as focal parameters, then replacing p by αp (while keeping U fixed) would no longer rescale the right-hand side of Eq. (28) by α, since the advection term would remain unchanged while the diffusion and source terms would change by a factor α. Intuitively, this would correspond to a rescaling of time (change in pace) for only a subset of processes (diffusion and reactions), while keeping time unscaled for advection, thus changing the balance between these three processes and ultimately altering the steady state.

The boundary concentrations A_j and B_j cannot be included as focal parameters, since this would violate the invariance condition for steady states. Indeed, changes in A_j or B_j would inevitably change the system’s steady state, at least at the boundary and in its vicinity, regardless of how the other model parameters are varied. For the stoichiometric coefficients ({mathbb{S}}_{{{j}}{{r}}}) the situation is more subtle. For example, if one or more stoichiometric coefficients were to be included in p along with the D_j, v_j, β_r, then replacing p by αp while keeping U fixed would generally not rescale the right-hand side of Eq. (28) by that same factor. However, if one were to omit all affinities β_r and instead include all stoichiometric coefficients in p, while also keeping the diffusivities D_j and advection speeds v_j, then one would again obtain a set of focal parameters. This is intuitive and may even be considered trivial, since uniformly rescaling all stoichiometric coefficients is equivalent to uniformly rescaling all reaction affinities, but it serves to highlight the fact that a model may in principle exhibit more than one possible set of focal parameters. A similar ambiguity exists in classical MCA, where one could in principle focus on stoichiometric coefficients instead of enzyme concentrations. In practice this is virtually never done, since stoichiometric coefficients are generally considered fixed and/or normalized according to some convention, for example such that electron acceptors have a coefficient of 1, or such that the coefficients of a reaction are co-prime.

Following the earlier discussion on grouping parameters, for simplicity I henceforth group the diffusivities of each metabolite across space. Specifically, I consider the average diffusivities ({tilde{D}}_{j}:= frac{1}{b-a}int_{a}^{b}{D}_{j}(z),dz), while keeping the location-dependent normalized diffusivities ({d}_{j}(z):= {D}_{j}(z)/{tilde{D}}_{j}) fixed. Due to the strong relationships between advection speeds, we group all advection speeds across space and across all metabolites into a single average, (tilde{{{{rm{v}}}}}), while keeping the metabolite- and location-dependent normalized advection speeds ({v}_{j}(z):= {{{{rm{v}}}}}_{j}(z)/tilde{{{{rm{v}}}}}) fixed. Thus, the vector of focal parameters considered for MCA is:

Control coefficients: We are now ready to examine the FCCs of the mean rate of methane oxidation and the CCCs of the mean concentration of hydrogen sulfide with respect to the above focal parameters (Eqs. (15) and (10), respectively). Numerically, the control coefficients can be estimated via centered difference quotients, i.e., by varying one parameter at a time by a small amount and examining the corresponding small changes in the steady state (details in Supplementary Discussion 2). FCCs for the mean methane oxidation rate (Fig. 1M) are all positive, and — due to the FCC summation theorem — can thus be interpreted as fractions of overall control. By far the most important FCCs are those associated with methane diffusivity and sulfate diffusivity, which together contribute about 99.8% of the control, while advection speed contributes about 0.2% of the control. In contrast, the FCC associated with the reaction’s affinity (β_AOM) is essentially zero (up to numerical error). Thus, the rate of methane oxidation appears to be nearly entirely controlled by the molecular diffusion rates of methane and sulfate, with downward advection contributing a detectable but minuscule part, and reaction affinity exerting practically no control at all. Ultimately, this dominance of transport as a controlling factor stems from the vast difference in time scales involved in vertical transport across the sediment column and the chemical reactions. Under these conditions, AOM is restricted to a narrow depth zone, whose location and total reaction rate is largely determined by the reaction stoichiometry, the imposed boundary conditions and the diffusivities of methane and sulfide^8,9. Hypothetically, if reaction rates were much slower or diffusion much faster, to the point where their time scales become comparable, reactions would not be confined to a narrow zone and their rates would be partly controlled by kinetic parameters. Indeed, repeating the above analysis with 10⁵ times lower affinities, substantially widens the zone where each reaction occurs, and makes FCCs associated with affinities comparable to those associated with diffusivities (Fig. S1).

CCCs for mean H₂S concentration (Fig. 1N) are largely dominated by the diffusivities of methane and sulfate (positive control) and the diffusivity of H₂S (negative control), followed by the much weaker (but detectable) controls associated with the diffusivity of Fe²⁺ (negative) and advection speed (positive). CCCs associated with reaction affinities have negligible control on the mean H₂S concentration. Thus, the extent of H₂S accumulation is largely determined by a balance between methane and sulfate transport (which controls H₂S production rates) and H₂S transport away from the area of production. Note that while some CCCs are positive and some are negative, as predicted by the concentration control summation theorem positive and negative CCCs ultimately cancel each other out.

Example: Fixed nitrogen loss in a seasonal oxygen minimum zone

As a second detailed example, I examine a 1-dimensional PDE model for microbial populations and microbially driven redox reactions in a seasonal oxygen minimum zone (OMZ) in Saanich Inlet, off the coast of Vancouver Island^23,24. OMZs are regions in the ocean water column where oxygen-poor conditions favor microbial anaerobic metabolisms, which commonly result in the conversion of bioavailable (fixed) nitrogen to dinitrogen (N₂)³⁸. OMZs thus act as sinks for bioavailable nitrogen in the ocean, with major implications for marine ecosystems³⁹. The model considered here was previously described and fitted to data by Louca et al.²², and hence here only an overview is given. The model is constructed to reproduce the quasi-steady state established in the Saanich Inlet water column in early 2010, during a period of intense stratification. The model tracks the concentrations of dissolved ammonium (NH4⁺), oxygen (O₂), nitrate (({{{{rm{NO}}}}}_{3}^{-})), nitrite (({{{{rm{NO}}}}}_{2}^{-})), hydrogen sulfide (H₂S), nitrous oxide (N₂O), sulfate (({{{{rm{SO}}}}}_{4}^{2-})) and dinitrogen (N₂), between depths a = 100 m and b = 200 m. The dynamics of these metabolites are coupled through 5 microbially driven reactions: Aerobic remineralization of particulate organic matter (AEROM):

aerobic ammonium oxidation to nitrite (AAO):

aerobic nitrite oxidation to nitrate (ANO):

anaerobic ammonium oxidation via reduction of nitrite (known as anammox⁴⁰):

partial denitrification to nitrite and nitrous soxide via H₂S oxidation (PDNO):

and reduction of nitrous oxide to N₂ coupled with H₂S oxidation (RNO):

The concentration profile of organic matter is explicitly fixed in the model. The above reactions are largely fueled by H₂S and ({{{{rm{NH}}}}}_{4}^{+}) diffusing upward from the sediments, and O₂ and sinking organic matter from shallower waters. All metabolites are subject to vertical eddy diffusion, which varies across depth depending on the local buoyancy frequency⁴¹. At steady state, each metabolite’s concentration thus satisfies the following PDE:

where as before z denotes depth, D is the depth-dependent eddy diffusivity (m² ⋅ day⁻¹), ({mathbb{S}}) is the stoichiometric matrix of the reactions and R_r is the rate of the r-th reaction (mol ⋅ L⁻¹ ⋅ d⁻¹). In contrast to the earlier SMTZ example, this model explicitly considers the population dynamics of microbial species catalyzing the reactions, with microbial growth being driven by the energy released through the reactions. Each reaction rate R_r is thus proportional to the cell density of a unique microbial species catalyzing the reaction, henceforth denoted W_r and measured in cells per L. Per-cell reaction rates are assumed to follow first-order (mass action) or Monod kinetics in each substrate, with some reactions further inhibited by some metabolites:

where L_r,1 and L_r,2 specify the substrates with first order and Monod order kinetics, respectively, for the r-th reaction, I_r specifies the inhibitors of the reaction, K_r,j is the half-saturation constant for substrate j in reaction r, ({K}_{r,j}^{star }) is the half-inhibition constant for inhibitor j in reaction r, β_r is the cell-specific affinity of the r-th reaction, and F_r is a thermodynamic potential factor that depends on metabolite concentrations and lowers the reaction rate when the latter is close to thermodynamic equilibrium⁴² Eq. 16. Note that here the affinity β_r is defined as the mass-action kinetic constant in the limit where all substrate and inhibitor concentrations approach zero. Cell death is modeled as an exponential decay process, and cells are subject to vertical eddy diffusion and sinking⁴³. Thus, at steady state each cell density W_r satisfies the PDE:

where δ_r is the constant exponential death rate for the r-th species (in units 1/day), Y_r is the biomass yield for that reaction (dry biomass produced per mole reaction flux), M is the average dry cell mass and v_r is the cell sinking speed. The yield Y_r is not a fixed parameter, but an empirically fitted linear function of the reaction’s Gibbs free energy⁴⁴, and thus depends on metabolite concentrations and the reaction stoichiometry. Note that the model’s dynamic variables include both U and W, that is, the system’s state is defined by the tuple (U, W).

Boundary conditions for metabolite concentrations at the top are all of fixed-value type, i.e., U_j(a) = A_j, where A_j is a constant. At the bottom, some metabolites have fixed-value conditions while others have zero flux conditions:

For cell densities, boundary conditions are either fixed values, or expressed in terms of fixed relative slopes, i.e.:

Additional details, including parameter values, boundary conditions and the eddy diffusivity profile, can be found in the Supplement of²². An overview of the steady state predicted by the model for default parameter values, including metabolite concentration and cell density profiles across depth, is shown in Fig. 2. While microbial cells responsible for aerobic reactions (AEROM, AAO and ANO) are mostly abundant toward the upper end of the depth interval, anaerobic cells engaged in nitrogen reduction thrive at greater depths, causing the formation of a sharp sulfide-nitrate interface at around 130 m. For an extensive discussion of this steady state see ref. ²².

Focal parameters: Similarly to the earlier SMTZ example, let us focus on the model’s rate-like parameters, i.e., the diffusivity D, the cell sinking speeds v_r, the reaction affinities β_r and the per-capita cell death rates δ_r. A uniform rescaling of these parameters by a constant factor (while keeping U and W fixed) rescales the right hand side of the PDEs (Equations (36) and (38)) by that same factor, it rescales the zero-flux boundary conditions (Eq. (39)) by that factor, and leaves the fixed-value and fixed-relative slope boundary conditions (Eq. (40)) unchanged. Thus, this set of parameters satisfies the homogeneity condition for F_ω in Eq. (9) and consequently the invariance condition for steady states (Eq. (7)). Uniformly rescaling these parameters (while keeping U and W fixed) also rescales each reaction rate R_r by the same factor, and thus reaction rates (and their averages across space) satisfy the homogeneity condition for fluxes in Eq. (13). A similar conclusion holds for the production and consumption rates of metabolites, as well as for microbial cell generation and death rates. Following the earlier discussion on grouping parameters, for MCA it is simpler to consider the average diffusivity (tilde{D}:= frac{1}{b-a}int_{a}^{b}D(z),dz), while keeping the location-dependent normalized diffusivity (d(z):= D(z)/tilde{D}) fixed. Thus, the vector of focal parameters is:

Control coefficients: While the model makes predictions about many different fluxes and coupled reactions, for MCA I henceforth focus on those directly involved in fixed nitrogen loss. Specifically, I consider the FCCs of the mean anammox, PDNO, and RNO reaction rates, the FCCs of the mean rate of N₂ production, the CCCs of the mean cell densities associated with anammox, PDNO, and RNO (Fig. 2M, N), as well as fluxes of nitrogen compounds across the top and bottom boundaries (100 m and 200 m depth). In contrast to the earlier SMTZ example, a small number of FCCs is actually negative, indicating that some focal parameters have a negative impact on process rates. The most negative FCC (-0.069) was obtained for anammox with respect to the per-capita death rate of AEROM cells. A likely explanation is that a higher mortality of AEROM cells reduces AEROM cell densities and thus AEROM reaction rates, in turn increasing the amount of oxygen that is available for aerobic ammonia oxidation and aerobic nitrite oxidation, both of which compete for nitrogen substrates with anammox, thus in turn reducing anammox rates. That said, negative FCCs are all relatively weak (i.e., close to zero), and thus — due to the FCC summation theorem — one can still approximately interpret the remaining non-negative FCCs as fractions of overall control. By far the strongest FCCs are those associated with the mean diffusivity ((tilde{D})), attaining values close to 1 in all cases. This suggests that eddy diffusion (turbulent mixing) was by far the most important modeled rate-limiting process for fixed nitrogen loss in Saanich Inlet during the considered time period. In contrast, FCCs associated with reaction affinities (β_r) are all close to zero, which implies that changes in the considered reaction affinities would only have a minimal impact on fixed nitrogen loss. Similar results are also obtained when focusing on the mean reaction rates around the redox transition zone (120–140 m, Fig. S2A). This suggests that for purposes of predicting overall fixed nitrogen loss rates in oxygen minimum zones it is more important to accurately characterize hydrodynamic mixing processes than the kinetics of supply-limited reactions. It should be noted that in this model the cycling of organic matter is not modeled fully mechanistically, since the depth-profile of organic matter was fixed and the production kinetics of organic matter are not considered at all. Organic matter remineralization, especially extracellular hydrolysis of organic macromolecules, is often considered the rate-limiting step in microbial communities^45,46, and in the present analysis the affinity and death rate of AEROM did have the strongest control on N₂ production among all reactions (Fig. 2). Further investigations are thus necessary to elucidate the control exerted by carbon cycling processes on fixed nitrogen loss.

Interestingly, the distribution of FCCs across parameters becomes more complex when one considers fluxes of nitrogen compounds (N₂, N₂O, ({{{{rm{NO}}}}}_{2}^{-}), ({{{{rm{NO}}}}}_{3}^{-}) and ({{{{rm{NH}}}}}_{4}^{+})) across individual boundaries (Fig. S2C). While the flux of ({{{{rm{NH}}}}}_{4}^{+}) across the bottom boundary is largely controlled by diffusivity (FCC 0.998, all other FCCs close to zero), for all other examined compounds FCCs cover a variety of values, including values far from 0 and 1 and even strongly negative values. For example, for N₂ flux across the top boundary, FCCs range from −0.65 and −0.58 with respect to the affinities of anammox and ANO, respectively, to 0.58 with respect to the affinity of AAO and up to 1.14 with respect to the diffusivity. A similarly complex partitioning is seen for N₂ flux across the bottom boundary. This suggests that the partitioning of N₂ release across the top and bottom boundary depends in complex ways on both transport and reaction kinetics, despite the fact that their sum (equal to total N₂ production) is almost entirely controlled by diffusivity alone (Fig. 2M).

CCCs for mean cell densities display a somewhat more diversified dependency (Fig. 2N). While CCCs associated with diffusivity (leftmost column in Fig. 2N) are still the strongest by magnitude, substantial control is also exerted by per-capita death rates (δ_r) and to a lesser extent by cell sinking speeds (v_r). On the other hand, reaction affinities exert almost no control on cell densities. These observations may come as no surprise, since steady state cell densities results from a balance between processes that add cells (mainly chemical reactions, whose rates are insensitive to affinities) and processes that remove cells (i.e., sinking and death).

Comparison to other common sensitivity analyses

It is worth noting where MCA stands in relation to other types of sensitivity analysis. Since FCCs and CCCs are based on partial derivatives of model responses with respect to focal parameters, MCA constitutes a local sensitivity analysis, that is, it only informs about the sensitivity of the system to small parameter changes. This contrasts MCA with global sensitivity analyses, for example those that “sweep” a wide range of parameter values⁴⁷ and Monte Carlo methods⁴⁸ (for a detailed comparison between local and global sensitivity analyses see refs. ^49,50). In fact, FCCs and CCCs are a special case of normalized local sensitivity coefficients^22,27, focused on specific types of observables (fluxes and concentrations) and focal parameters that control the system’s pace. The deliberate examination of a subset of pace-setting focal parameters, and the recognition of system responses with particular sensitivity properties (here, termed “fluxes” and “concentrations”), gives MCA additional structure and interpretability, including summation theorems that relate the sensitivity coefficients to each other, and the interpretation of FCCs as fractional contributions to flux control.

Classical MCA has been an invaluable tool in understanding the controls of cellular biochemical reaction networks modeled via ordinary differential equations. The generalization presented here covers a much wider range of reactive systems, such as microbial communities coupled to redox chemistry in water columns, sediments or bioreactors, elemental cycles at regional or planetary scales, and even ecological food webs (Supplementary Discussion 1). Mathematically, this generalized MCA can accommodate a broad range of model structures, including nearly arbitrary combinations of ordinary differential, partial differential, integral, and algebraic equations. Indeed, the two critical assumptions made are the invariance of steady states (Eq. (7)) and the 1-homogeneity of fluxes (Eq. (13)) to a uniform rescaling of the focal parameters, regardless of the precise model structure.

The focal parameters defined here present an intuitive generalization of the enzyme activities considered in classical MCA. As illustrated through multiple examples, focal parameters will often correspond to the rate-like parameters (and inverses of time-like parameters) in a model, which define the time scales of various processes and the “pace” of a system’s dynamics. In that sense, control coefficients formalize the dependency of a system’s emergent properties on the time scales of individual constituent processes, such as chemical reactions, microbial population growth or death, diffusive and advective transport, outgassing at water-air interfaces, or hydraulic renewal in bioreactors. The concentration and flux control summation theorems, in turn, reveal fundamental relationships between control coefficients; for example, since all CCCs must sum to zero, negative CCCs inevitably imply the existence of positive CCCs as well, and increasing some CCCs will necessarily decrease other CCCs.

In both of the detailed simulation examples discussed, the sulfate-methane transition zone in Black Sea sediments and a seasonal oxygen minimum zone in Saanich Inlet at quasi-steady state, physical transport processes across space were found to exert by far the greatest control on the fluxes of interest. This suggests that greater emphasis should be put on characterizing the transport processes in such systems, compared to the chemical and microbial kinetics. These examples demonstrate the great potential that MCA has for determining the rate-limiting processes in natural and engineered reactive systems.

While I focused on systems at steady state, it is expected that MCA as presented here is also approximately applicable to systems at quasi-steady state. Here, quasi-steady state is taken to mean that the system quickly converges to and closely follows a state (overline{{{{bf{U}}}}}) that is solely determined by the model parameters, with the latter either being constant or changing over much longer time scales than the time it takes for the system to converge to (overline{{{{bf{U}}}}}). In that case, predicted concentrations and fluxes, as well as their control coefficients, would be time-dependent quantities determined by the current parameter values. How far from steady state the theory presented here starts to fail remains an open and complex question. It might also be possible to rigorously generalize MCA to systems far from steady state, albeit at the risk of complicating interpretations. Indeed, classical MCA has been generalized to time-dependent biochemical networks^30,32, but its adoption has been limited compared to steady state MCA, partly due to the added complexity.

The predictive power of classical MCA and its potential to yield practical insight into cell metabolism have been confirmed through numerous experiments and successful bioengineering outcomes guided by MCA^18,19. By analogy, I anticipate that the potential of the generalized MCA presented here will become further apparent through successful predictions of geochemical/ecological flux changes in response to specific environmental shifts, as well as future eco- and geo-engineering successes guided by MCA.

Additional details on the structure, parameterization and simulation of the Black Sea SMTZ model are provided in Supplementary Discussion 2. The Saanich Inlet oxygen minimum zone model is described in detail in ref. ²².