Organization of biogeochemical nitrogen pathways with switch-like adjustment in fluctuating soil redox conditions

Sanjay Lamba1, Soumen Bera1, Mubasher Rashid1, Alexander B. Medvinsky2, Gui-Quan Sun3, Claudia Acquisti4, Amit Chakraborty1,5 and Bai-Lian Li5 1School of Mathematics, Statistics and Computational Sciences, Central University of Rajasthan, Bandarsindri, Ajmer, India 2Institute of Theoretical and Experimental Biophysics, Pushchino 142290, Russia 3Department of Mathematics, Shanxi University, Taiyuan, People’s Republic of China 4Institute for Evolution and Biodiversity, WWUMuenster, Germany 5Ecological Complexity and Modeling Laboratory, University of California, Riverside, CA 92521-0124, USA


Introduction
This Auxiliary Material includes five supplementary tables (Table S1-S5), fifteen supplementary figures (Figs. S1-S15) and ten supplementary texts(Text S1-S10). Detailed methods and mathematical techniques are described in the supplementary texts; all the simulated model results in the form of tables and figures are presented in the auxiliary figures and tables. This material presents how the nitrogen biochemical system runs and is regulated over time and how it is responded to external ammonium supply.

S1 Nitrogen Biochemical Network
The nitrogen biochemical network is constructed with the use of Kyoto Encyclopaedia of Genes and Genome database (KEGG,www.genome.jp/kegg). It consists of twelve commonly occurred biochemical pathways that serve to biochemically process nitrogenous metabolites and transfer the product into next reaction of the network ( Figure S1 ).
The network nodes are representing nitrogenous metabolites within a closed circle and a directed edge between the nodes symbolize enzymatic reaction that biochemically processes the metabolites and extract nitrogenous form different from the substrate. A rectangular box on the edge denotes biological structure excreted by a variety of soil microbes. It can be noticed that there are multiple biochemical pathways (edges) that catalyze by different biological structure between the same substrate and product of microbial actions, forming a complex biochemical nitrogen network. In the network, ammonium and nitrite act as network hubs connecting relatively a large number of nodes in the network. The graph representation of this network asserts that it forms a connected graph, as there exists a directed path between any pair of nodes. However, species graph representation of the associated biochemical reactions ensures that the nitrogen biochemical network cannot form a monotone dynamical system because of multiple biochemical pathways between a pair of nodes.
Table S1 is list the biochemical reaction pathways symbolized by r and the associated biological structure catalysing the pathway reaction.

S2 Mass Balance Equations of the Nitrogen (N ) Biochemical Network
The mass balance equations are formulated for the nitrogen biochemical reactions network which describe the dynamics of all nitrogenous metabolite concentrations and illustrate the steady-state relationship among the associated pathways. The equation is given by, where d X d t denotes time derivative of the metabolite concentration, X in the network, describing instantaneous change in simultaneously occurring biochemical reactions. S is the stoichiometric matrix representing the network topological structure, in which metabolites are represented by the nodes and the directed edges symbolize unidirectional biochemical reactions. The substratedependent enzymatic reaction rate, v (X ) , is described by the irreversible Michaelis-Menten kinetics. The stoichiometric matrix S is given by S = r 1 r 2 r 3 r 4 r 5 r 6 r 7 r 8 r 9 r 10 r 11 r 12  A system of ordinary differential equations describing time-dependent simultaneous changes of metabolite concentration x 1 (t ), x 2 (t ), x 3 (t ), x 4 (t ), x 5 (t ), x 6 (t ) and x 7 (t ) can be constructed from the vector S.v as follows: x 3 x 4 x 5 x 6 Hydroxylamine N H 2 OH denote the molar concentrations of metabolites. And the constants k 1 , k 2 , k 3 , k 4 , k 5 , k 6 , k 7 , k 8 , k 9 , k 10 , k 11 and k 12 are the reaction rate coefficients for the biochemical reactions r 1 , r 2 , r 3 , r 4 , r 5 , r 6 , r 7 , r 8 , r 9 , r 10 , r 11 and r 12 respectively (Tables S1 & S2). As we are interested in simulating the qualitative dynamics of the network, we assume, for simplicity, that the associated reaction rates, v i (X ) , has irreversible Michaelis-Menten form, k i X k m i +X (i = 1, 2, . . . , 12), with the half-saturation co-efficient k m i . Therefore, the transient behaviour of the metabolite concentrations are described by the following system of ordinary differential equations: Numerical simulation of the mass balance equations illustrates that the system can eventually reach to a steady-state at which each reactant maintains a fixed concentration determined by the reaction constants and the concentration of other metabolites ( Figure S2). In particular, we observe that the transient dynamics of ammonium and nitrate concentration follow very different trajectories. (Figure S3), While the ammonium concentration increases (decreases), the nitrate concentration decreases (increases) and then both concentrations finally attain its steady-state at which ammonium concentration remains higher (lower) than the nitrate concentration. On the other hand, while both the concentrations initially decrease by maintaining higher concentration of ammonium, the nitrate concentration jumps to an elevated level after a certain period of time and then concentration difference shrink gradually over a longer period of time. This qualitative behaviour is, therefore, indicating that nitrogen biochemical network involves a switching mechanism that flips the steady-state levels of inorganic ammonium and nitrate-the most common bioavailable forms of nitrogen.

S3 Response of N −Biochemical Systems to N −Fertilization
We have simulated this effect by adding ammonium as pulse input. Similar in-silico experiment can be carried out for the nitrate as well. To correctly locate the departure point from the biochemical steady-state, we have designed two different in silico experiments: (a) Ammonium is supplied at very beginning of the evolution of N −biochemical system, and (b) Ammonium is supplied at a transient state away from the steady-state. For understanding the qualitative characteristics of the departure point, four dynamic regimes are defined following the flipping behavior mentioned in the previous section:-I: Ammonium concentration is more than the 5 times of nitrate concentration (i.e., ammonia rich regimes), II: Low availability of both ammonium and nitrate (i.e., nitrogen poor regimes), III: Nearly identical concentration of ammonium and nitrate, and IV: Nitrate concentration is more than the 5 times of ammonium concentration (i.e., nitrate rich regimes). We have then changed the levels of pulse input of ammonium from near-to-zero to 15 fold increased level (Table S3, Figures S4 & S5).

S4 Design of Monotone Dynamical Systems
Smith [1], Hirsch and Smith [2] have introduced the idea of monotone dynamical systems with the properties of robust dynamical stability and predictability of responses to perturbations. Basic qualitative feature of this system is that its dynamical behaviour depends only on associated network structure not upon the precise values of the parameters or even the forms of reactions. This system describes the evolution of states which are time-dependent vector where component x i represents concentration of chemical species such as metabolities. Such reaction system has often been represented as autonomous system of ordinary differential equations with the form: where v (X ) is a m−dimensional vector of reactions and S is an n × m matrix called stoichiometric matrix.
A graph based approach has been developed to construct monotone dynamical system associated with a biochemical network. Species graph G represents the topological structure of the biochemical network, with all the nodes correspond to chemical species and the directed edges between the nodes symbolize the biochemical reaction. A dynamical system is then said to be monotone if there exists at least one consistent spin assignment for its associated species graph G. Continuous-time monotone systems have convergent behaviour; they cannot admit any possible stable oscillations. If there is only one steady-state and solutions are bounded, the every solution converges to this unique steady-state. Instead, if there are multiple steady-states, the Hirsch Generic convergence theorem asserts that generic bounded solutions of strongly montone system must converge to the set of steady-states. In particular, no "chaotic" or any other "strange" dynamics can occur (see the details in Sontag [3]).
Positive feedback interconnections between monotone systems with the output of one enter as input to other establish bistability. Besides having intuitive proof of bistability for simple production/degradation system with sigmoidal characteristics, this phenomenon has well demonstrated for Mitogen-activated protein kinase (MAPK) cascades [4].
We have adopted Montone dynamical system for decomposing the biochemical nitrogen network. A spin assignment and consistency based graphical approach has been implemented to derive a monotone subsystem, which says that a dynamical system is monotone if and only if there exists at least one consistent spin assignment for its associated reaction graph. In addition, there will be a consistent spin assignment if and only if every undirected loop in the graph has an even number or zero of negative sign. For constructing a monotone dynamical system ( Figure S6), following rules are considered that are associated with a monotone graph: where the positive terms in the equations describe the formation and the negative terms describe the degradation rate of metabolites. At a steady-state, degradation rate equates with the formation rate which is given by It shows that steady states lies on the curve ux = (k m 1 )(k m 2 ).( Figure S7) Bistability is often detected in biochemical systems that contain a positive feedback loop or double-negative feedback loop. A bistable system exhibits flipping behaviour between two steady states while passing through an unstable steady-state. It has been observed among various biochemical systems with different orders that the presence of at least one positive feedback loop is necessary for the emergence of bistability (Thomas [6]). However, mere presence of positive feedback loop does not guarantee bistability. Indeed, this is one of the common technical problems for detecting bistability in higher order system where phase plane analysis often fails. Using the classic theory of monotone dynamical systems, Angeli [4] has established an elegant method for the analysis of bistability in positive feedback systems of arbitrary order. It showed that if the openloop, feedback-blocked system is monotone and possesses a sigmoidal response characteristic, then the feedback system always produces bistability for some range of feedback strength.
We have employed this method for explaining observed bistability in ammonium formation in the N −biochemical system. Using the graphical method (Sontag [3]), we have designed two monotone subsystems, termed as ammonium-source and sink systems, embedded in the N − biochemical system. Each subsystem produces monostable ammonium state. The source system receives nitrate as input and produces ammonium as an output, whereas the sink system acts in reverse way. A positive-feedback connection between these two subsystems is given by the ammonium oxidation pathways which is mediated by amo and hao encoded enzymes converting ammonium to hydroxylamine and then to nitrite . Switching from the lower to higher ammonium state is accompanied by transiently augmenting ammonia oxidation that rapidly reduces ammonium concentration. Transient reduction of inorganic ammonium concentration promotes activities of assimilatory nitrate reductase by relaxing ammonium's inhibitory effects. Once it crosses the inhibitory threshold, assimilatory pathways contribute more to total ammonium formation relative to DNRA.

S5 Ammonium Source Systems
Associated stoichiometric matrix, S sour ce , describing the N H + 4 −source topological structure (Figure S8) is given by Montone dynamical system for N H + 4 −source systems is, therefore, written as:

S6 Ammonium Sink Systems
Associated stoichiometric matrix, S si nk , describing the N H + 4 −sink topological structure ( Figure  S11) is given by Monotone dynamical system for the N H + 4 −sink pathways is, therefore, written as: x 6 k m 9 + x 6 − k 10 x 7 k m 10 + x 7 . 6 S7 Simulating Responses of N H + 4 −Source and N H +

−Sink Systems to External Step Inputs
We have simulated the responses of N H + 4 −source and N H + 4 −sink systems to step inputs in the formation of ammonium, with the underlying assumption that both the system will receive ammonium from external sources either through N fertilization or biological N fixation that effect microbial nitrogen pathways. Tables (S4 & S5) and the associated diagrams shows the responses of N H + 4 source-sink systems (Figures S14 & S15). Differential response patterns of these two systems indicate disparities in the roles of N H + 4 −sink and N H + 4 −source subsystems in the nitrogen biochemical network.

S8 Appendix A: Dependent and Independent Biochemical Reactions
With the use of standard matrix method, we have been determined the dependence and independence of biochemical reactions represented by the equation, where I is Identity Matrix.
We then apply elementary row operations both sides of the above equation until stoichiometric matrix S is reduced to echelon form R, where We have partition the M matrix in to U and V, along the same row line as the reduced echelon form.
x 3 x 4 x 5 x 6 x 3 x 4 x 5 x 6 Therefore the equation Multiplying out the lower partition, we obtain V d X d t = 0, ⇒ẋ 1 +ẋ 2 +ẋ 3 +ẋ 4 +ẋ 5 +ẋ 6 +ẋ 7 = 0, the independent equations are Since the bottom S 0 dependent rows can be derived by linear combinations of the top S R rows, we can define a matrix, the link matrix(L 0 ),that can carry out this operation: We know that where X i are the independent species and X d are the dependent species.
Equation (1) multiply by L 0 and subtract by equation (2) then we have Integrating above equation then we get x 3 x 4 x 5 x 6
Substituting these into the equations, we obtain following new differential equations in the basin of equilibrium point x * 1 , x * 2 , x * 6 .
Expanding R.H.S of equation (3) with the help of Taylor's Series about x * 1 , x * 2 , x * 6 and neglecting second and higher order terms, we geṫ Superscript * indicate values are at equilibrium point. Also Let us take From equations (4) and (5), we obtain Which is a system of linear equations. Therefore, for non trivial solution of the system, we have Simplifying above determinant, we get cubic equation in λ, which is known as characteristic equation of the system and is given as B * = −k 1 k 2 + (k 1 k 12 + k 11 k 12 + k 12 + k 2 ) x * 6 1 + x * 6 − (k 1 k 12 + k 11 k 12 + k 12 + k 2 ) 1 Now, solving (8) for λ, we get From equation (9), it is clear that all the eigenvalues are real and negative. Therefore, the equilibrium point is a stable node or sink i.e., asymptotically stable.

S10 Appendix C: Stability of N H + 4 −Sink Systems
With the assumption that intermediate N −forms hydroxylamine and nitrite transforms faster, the stability of the system is therefore determined by the dynamics of ammonium and nitrate only. Let X 1 and X 6 be the small perturbation to the steady state concentration x * 1 , x * 6 so that the co-ordinates of N H + 4 −sink metabolites and components of its rate of change are x 1 = x * 1 + X 1 , x 6 = x * 6 + X 6 andẋ 1 =Ẋ 1 ,ẋ 6 =Ẋ 6 respectively, where X 1 , X 6 andẊ 1 ,Ẋ 6 are very small quantities.
Substituting these into the equations, we obtain following new differential equations in the basin of equilibrium point x * 1 , x * 6 .Ẋ Expanding R.H.S of equation (10) with the help of Taylor's Series about x * 1 , x * 6 and neglecting second and higher order terms, we geṫ Superscript * indicate values are at equilibrium point. Also f 1 x * 1 , x * 6 = f 6 x * 1 , x * 6 . Let us take From equations (11) and (12), we obtain Which is a system of linear equations. Therefore, for non trivial solution of the system, we have Simplifying above determinant, we get cubic equation in λ, which is known as characteristic equation of the system and is given as Where Now, solving (15) for λ, we get From equation (16), it is clear that all the eigenvalues are real and negative. Therefore, the equilibrium point is a stable node or sink i.e., asymptotically stable.          Step Input (Red) Output No.
Upper Level Lower Level x 1 (Black) x 6 (Cyan)