Analytical singular value decomposition for a class of stoichiometry matrices

We present the analytical singular value decomposition of the stoichiometry matrix for a spatially discrete reaction-diffusion system on a one dimensional domain. The domain has two subregions which share a single common boundary. Each of the subregions is further partitioned into a finite number of compartments. Chemical reactions can occur within a compartment, whereas diffusion is represented as movement between adjacent compartments. Inspired by biology, we study both 1) the case where the reactions on each side of the boundary are different and only certain species diffuse across the boundary as well as 2) the case with spatially homogenous reactions and diffusion. We write the stoichiometry matrix for these two classes of systems using a Kronecker product formulation. For the first scenario, we apply linear perturbation theory to derive an approximate singular value decomposition in the limit as diffusion becomes much faster than reactions. For the second scenario, we derive an exact analytical singular value decomposition for all relative diffusion and reaction time scales. By writing the stoichiometry matrix using Kronecker products, we show that the singular vectors and values can also be written concisely using Kronecker products. Ultimately, we find that the singular value decomposition of the reaction-diffusion stoichiometry matrix depends on the singular value decompositions of smaller matrices. These smaller matrices represent modified versions of the reaction-only stoichiometry matrices and the analytically known diffusion-only stoichiometry matrix. Our results provide a mathematical framework that can be used to study complex biochemical systems with metabolic compartments. MATLAB code for calculating the SVD equations is available at \url{www.github.com/MathBioCU/ReacDiffStoicSVD}.


Introduction.
In stoichiometric network analysis the mass balance equation for a reaction-only system is written as follows (1.1) dw dt = S r f where w is a species concentration vector, S r is the stoichiometry matrix, and f is a vector of reaction fluxes [2].We use the subscript r to refer to a stoichiometry matrix that only describes reactive processes.Although the flux vector f is a function of the species concentration, the formulation given by (1.1) avoids assumptions about the form of the kinetic equations that relate the fluxes to the species concentration (e.g., mass-action [14] or Michaelis-Menton kinetics [5]).The stoichiometry matrix contains information about the species involved in each reaction.As a simple example, consider the following set of reactions: Here, species A is produced, transitions into species B, and species B decays.The stoichiometry matrix for this example system is The first row of S r corresponds to species A and the second row corresponds to species B. Each of the three columns correspond to the three reactions, respectively.
The analysis of S r provides information on structural properties of the system without requiring kinetic information.In particular, the singular value decomposition (SVD) of S r provides information on systemic properties, including decoupled eigenreactions (i.e., linear combinations of species that are moved by linear combinations of fluxes), conservation relations, and fluxes that can exist in the system under steady-state conditions [11].This type of analysis can be used to determine hidden relationships in a network and compare biochemical properties amongst different organisms [4,10].
Here, our goal is to derive the SVD of a stoichiometry matrix that, in addition to the reactions, includes information on the spatial properties of a system.Specifically, we define the stoichiometry matrix for a one dimensional spatially discrete system by considering both the reactions in each spatial compartment as well as the movement of species between adjacent compartments.We refer to this matrix as the reaction-diffusion (RD) stoichiometry matrix and write it using the reaction-only and diffusion-only stoichiometry matrices, i.e. S r and S d , respectively.The reaction only stoichiometry matrix is as described in (1.3) and the diffusion only stoichiometry matrix can be thought of as representing a single species diffusing through space.As an example, in a system with one diffusing species, three spatial compartments, and homogeneous Neumann boundary conditions, the diffusion-only stoichiometry matrix is (1.4) Analagously to S r , each row in S d corresponds to the species in each of the three compartments.The first and last column of S d correspond to species movement across the boundary of the domain and, for this example, contain only zeros due to the homogenous Neumann boundary conditions.The middle columns represent the movement of the species between adjacent compartments.Using the reaction-only and diffusion-only stoichiometry matrix definitions, in a system with both reactions and diffusion where there are n spatial compartments and m species that freely diffuse through space, the stoichiometry matrix S is where ⊗ represents the Kronecker product [8], γ > 0 describes the relative rate of reactions to diffusion, and I a is the identity matrix of size a.Here, the S r ⊗ I n block represents the reactions occurring in each compartment, whereas the I m ⊗ S d block represents diffusive movement.The Kronecker product has previously been used to compactly represent diffusion for the spatially discrete reaction-diffusion ODE system [1,3].We previously developed criteria to guarantee a version of this ODE system is bounded for all time [15].Here, we instead use (1.5) to study the spatially discrete system in the context of stoichiometric network analysis.We write both the reactive and diffusive terms using a Kronecker product formulation as this will simplify the SVD derivation.
In this paper we will consider a more general form of (1.5) where, in addition to diffusion, there is a spatial barrier in the system that divides the 1D domain into two subregions.We consider this class of systems because it allows our results to be applied to study, for example, the effect of metabolic compartmentalization within a cell.We will use concepts from linear perturbation theory [6] to derive the approximate SVD in the limit as diffusion becomes much faster than reactions.We additionally consider the special case where diffusion of all species is allowed freely throughout the domain, i.e., where the stoichiometry matrix can be written as given by (1.5).We show that, for this scenario the SVD becomes exact for all values of γ.The derived SVDs for the system with and without a spatial barrier depend on the SVDs of smaller matrices, such as the reaction-only stoichiometry matrix.
To help provide structure and guide our argument, in Section 2 we chose to present the main result first (see Theorems 2.1 and 2.2).We then provide a more complete set of definitions and notation in Section 3.This includes a complete description of the system as well as definitions of matrices whose SVDs are used to write main result.In Section 4 we provide preliminary results that will be helpful for proving Theorem 2.1.In Section 5 we provide the complete proofs of Theorem 2.1 and Theorem 2.2.Finally, in Section 6 we provide some intuition for the SVD equations and discuss potential applications of this work.
2. System description and statement of main result.Here we provide a brief description of the system and state the main result.For a thorough description of the notation and definitions used see Section 3.
We consider a one dimensional, spatially-discrete, reaction diffusion system that is divided into two subregions.A subset of the species is allowed to diffuse between the two subregions, and we allow for different sets of reactions to occur in each region.We will consider three boundary conditions: no input/output fluxes, input/output fluxes at one boundary point, and input/output fluxes at both boundary points.As an example, biologically this system description might represent a radially symmetric cell, where the two subregions are the cytoplasm and the nucleus.
The stoichiometry matrix for this class of systems can be written as (2.1) where the first column block represents reactive processes and the second represents diffusive processes.Here, S r1 and S r2 represent the reaction-only stoichiometry matrices for each of the two subregions, S d ⊗ D + describes the diffusion of species that move across the entire domain (i.e., species that can cross the barrier between the two subregions), and (S d − H) ⊗ D − describes the diffusion of species that stay within a single subregion.The parameter γ > 0 represents the relative rate of reactions compared with diffusion.
In this section we present the SVD of the stoichiometry matrix given by (2.1) in the limit as diffusion becomes much faster than reactions, i.e., as γ → 0. Briefly, the main result depends on the SVD of smaller reaction-only and diffusion-only systems.This includes matrices that only involve reactive processes, which will be written using variations of S r (e.g., S r1 , S r2 ), and matrices that only involve diffusive processes, which will be written using variations of S d (e.g., S d1 , S d2 ).
Our general notation for writing down the SVD of S• ∈ R s1×s2 will be as follows: (2.2) We will refer to the rank of S• as q• and the size of the nullspace as q•.In some cases the singular vectors will be divided into two components (e.g., ).With this SVD notation in mind, we next state the main result of the paper.Although the complete definitions and notations are not given until Section 3, it is possible to immediately see that the SVD depends only on SVDs of variations of stoichiometry matrices for the reaction-only and diffusion-only systems.
Theorem 2.1.As γ → 0 the unsorted SVD of S, as given by (2.1), is where the singular vectors that have nonzero singular values are given by six components, Ûi , Vi , such that A basis for the left nullspace of S is where and a basis for the (right) nullspace of S is Note that the horizontal dashed lines used in the definition of the right singular vectors separate the vectors into components that correspond to the reactive fluxes (above dashed line) and diffusive fluxes (below dashed line).The proof of this theorem is given in Section 5.
We have defined the SVD in Theorem 2.1 to be applicable for all three boundary conditions.Note that Û6 , V6 , Ȗ2 , and V3 are only nonempty for homogeneous Neumann boundary conditions (i.e., no input/output fluxes) and Û3 , V3 , Ȗ1 , and V1 are only nonempty when there is an input/output flux at a single boundary point.
The results given in Theorem 2.1 are simplified significantly when we consider systems that only have one region and spatially-homogeneous reactions.For such systems the stoichiometry matrix is simplified to and the SVD is given by the following theorem.
Theorem 2.2.The SVD of the stoichiometry matrix S, as given by (2.8) is 3. Notation and Definitions.Here we present notation and matrix definitions that are used to state and prove the main result.In Section 3.1, we provide basic notation for referring to matrices.In Section 3.2, we present definitions used to define the discrete reaction-diffusion system.In Section 3.3 we define sets of indices that will be used for defining the SVD.In Section 3.4, we define the set of stoichiometry-like matrices that are required for writing the SVD of the reaction-diffusion system.In Section 3.5, we provide notation that, in addition to (2.2), will be used to define the SVD of relevant matrices.Table 1 summarizes the notational defintions presented in this section.
3.1.Matrix notation.Matrices will be defined using uppercase letters (e.g., A) and sets of indices will be defined using calligraphic fonts (e.g., B).We will use A B to represent the columns of A whose indices are in the set B. If we refer to one column of a matrix (i.e., a column vector), we will typically use the lowercase letter and a superscript to refer to this column (i.e., the ith column of A will be written as a (i) ).One exception to these rules will be for any diagonal matrix of singular values Σ• and variations of this matrix.In this case, Σ B • will represent a square diagonal matrix containing the singular values whose indices are in B. Additionally (Σ•) n will represent the matrix Σ padded by zeros to make it size n × n.The ith diagonal element of Σ• will we written as σ (i) • .Throughout the paper, we will use I a to denote the identity matrix of size a × a.We will also use 0 to represent a matrix of zeros.For notational simplicity we omit the size of each zero matrix but note that it can be deduced from the notation.We will use ⊗ to represent the Kronecker product 1 and ⊕ to represent the Kronecker

Discrete reaction-diffusion systems.
We consider the discrete reactiondiffusion system on a one dimensional domain [0, n] that is partitioned into n equalsized spatial compartments.Let m denote the number of species (e.g., proteins or metabolites) in the system and p denote the number of reactions.We will allow for three different boundary conditions: homogeneous Neumann (no flux at both ends), Mixed (homogeneous Neumann at x = 0 and open at x = n), and Open (flux allowed at both ends).Note that both the reactive and diffusive fluxes can be either positive or negative.We define a positive diffusive flux as moving in the positive x direction.We will assume that all the species in the system diffuse at the same rate.
Within the domain there is a single barrier across which only a subset of species can diffuse.The barrier divides the system into two subregions where different reactions occur.Let S r1 ∈ R m×p1 and S r2 ∈ R m×p2 represent the stoichiometry matrices for the two subregions (i.e., p 1 reactions occur in the first subregion and p 2 reactions occur in the second).Note that the same reaction can occur in both regions.
We will let n 1 ∈ {1, ..., n − 1} denote the number of compartments in the first subregion and n 2 = n − n 1 denote the number of compartments in the second sub- 2 The Kronecker sum is given by where A is an a × a matrix and B is a b × b matrix.
region.Within this system, there are three diffusive processes: diffusion across the entire domain, within the first subregion, and within the second subregion.We define diffusion-only stoichiometry matrices for these three processes using S d ∈ R n×n+1 , S d1 ∈ R n1×n1+1 and S d2 ∈ R n2×n2+1 , respectively.For diffusion across the entire domain, we have that (3.1) where the values in the first and last column depend on the boundary conditions.
, we can relate S d with S d1 and S d2 as follows Next, we provide definitions used to identify the species that can and cannot diffuse between the two subregions.When defining parameters (e.g., sets, matrices), a subscripted + or − will imply a relationship with the set of species that can (+) or cannot (−) diffuse across the barrier.The set M + will contain indices for species that can diffuse across the barrier, whereas the set M − will contain indices for species that cannot diffuse across the barrier.Additionally, let m We can now write the equation for the spatially-discrete RD stoichiometry matrix, given by (2.1).For convenience we rewrite this equation below The parameter γ ≥ 0 represents the relative rate of the reactions compared with the rate of diffusion (i.e., if γ 1 the reactions are much faster than diffusion and, if γ 1, the reactions are much slower than diffusion).The first n 1 p 1 + n 2 p 2 columns of S correspond to the reactions occurring in each compartment.The final (n + 1)m columns correspond to the diffusion of species into or out of the domain as well as between adjacent compartments.

Additional spatially-dependent parameters.
Here we define the constants C 1 , C 2 , the set B, and sets denoted by variations of J .These parameters are only dependent on the spatial properties of the system (e.g., compartment number and boundary conditions), and are therefore unaffected if reactive properties (e.g., reaction number and stoichiometry) change.
The constants C 1 and C 2 depend on the boundary conditions and compartment numbers.We have that We will show in Lemma 4.1 that these constants relate the singular vectors for S d , S d1 and S d2 to one another.
Next, the set B is defined to contain indices that correspond to the columns of S d that are zero as well as the index of the column of S d that corresponds to the flux between the two subregions.Specifically, This set will be used to help define the nullspace of the RD stoichiometry matrix.Finally, we define the following index sets of singular values for the diffusion-only stoichiometry matrices (3.7) and the analogous index sets for only nonzero singular values We will use these sets to define how singular values repeat in the system when γ = 0. Understanding this property is a key step in proving Theorem 2.1.

Additional reaction-dependent stoichiometry-like matrices.
We refer to modified versions of the reaction-only stoichiometry matrices as stoichiometrylike matrices.In this section we will define the 4 + | Ĵ | stoichiometry-like matrices that are necessary for writing the SVD.These matrices are given as S r1,− , S r2,− , S r and S r+,j for j ∈ Ĵ .
The matrices S r1,− and S r2,− will represent subsetted versions of S r1 and S r2 , respectively, that only contain rows for species that cannot diffuse across the boundary.Next, we define a stoichiometry-like matrix that represents a merger of the two reaction-only stoichiometry matrices: .
To prove Theorem 2.1, we will need to consider the eigendecomposition of It can be shown that To obtain this equation, we use that I M+ m S r1,+ = D + S r1 and similar identities.Finally, for j ∈ Ĵ , define (3.14) S r+,j := |u Note that S r+,j for j ∈ Ĵ are the only stoichiometry-like matrices that depend on the spatial properties of the system.
3.5.Additional SVD notation for stoichiometry-like and the diffusiononly stoichiometry matrices.Generally, (2.2) will be used to write the SVDs of the stoichiometry and stoichiometry-like matrices.However, there are a few additional notational notes and one exception that will be discussed in this section.
First, the exception to this notational format will be for the left singular vectors of S r .Specifically, when considering the left nullspace of S r , we will exclude the space spanned by the following set of vectors We define Ȗr := span(null(S T , Ȗr,ex )) and U r := Ûr Ȗr .The reason for this will become clear in the proof to Theorem 2.1.
In some instances, we divide a given singular vector into two components.We will use a subscripted s 1 or s 2 to refer to portions of the singular vectors that correspond to processes that occur in the first or second subregion, respectively.Additionally, we wil use the subscript m 1 and m 2 to represent singular vectors that are divided into two subvectors of size m.More specifically, for the singular vectors of S d , we have that (3.16) where For the singular vectors of S r and the right singular vectors of S r+,j , we define where , and v r+,j,s2 ∈ R p2 .We will use the same notation to divide an entire set of right or left singular vectors into components.As an example, we have that , Ûd = Ûd,s1 Ûd,s2 , Ȗd = Ȗd,s1 Ȗd,s2 .
When considering the SVD of the diffusion-only stoichiometry matrices S d , S d1 and S d2 , the singular vectors and values can be written explicitly and depend on the specific boundary conditions (see Supplemental Material B).The rank of S d , given by q d , also depends on the the boundary conditions where and b 2 is as given in (3.1).This implies that the left nullspace, spanned by Ȗd , is empty for both Mixed and Open boundary conditions.
4. Preliminary Lemmas.In this section, we will provide preliminary lemmas that will be used to prove the main result.
First we consider the SVD of the diffusion-only stoichiometry matrices.In the following lemma we prove that, if a given singular value repeats across Σd , Σd1 , and Σd2 , then it must be in all of these matrices.That is, a singular value will occur in either one or all three matrices.Lemma 4.1.Consider a system with Zero Flux, Mixed, or Open boundary conditions and the singular values defined in the matrices Σd , Σd1 , and Σd2 .If a singular value is in two of these matrices then it is in all three.
For singular values that are in all three matrices, the corresponding singular vectors are related as follows: where j, j 1 , and j 2 are such that σ d1 .Additionally, the indices j, j 1 , and j 2 satisfy j 1 = C 2 1 j, j 2 = C 2 2 j, and j = j 1 + j 2 .The proof of this claim is given in Appendix D.
We next derive formulas for the dimensions of the four fundamental subspaces of S, as given by (2.1).This allows us to verify that the SVD has the correct number of singular vectors in each space.
Lemma 4.2.The rank of S, as given by (2.1), is Open.
The dimension of the nullspace is and the dimension of the left nullspace is Open.
Here we are using the rank and nullspace size of S r1,− and S r defined in (3.9) and (3.10).We omit the proof of this claim but note that it involves a sequence of row and column operations on S.

Singular value decomposition derivation.
In this section, we present the proofs for Theorem 2.1 and 2.2.Recall in Theorem 2.1 we provide the approximate SVD for a system with a barrier, whereas in Theorem 2.2 we consider a system without a barrier and derive an exact SVD for all relative diffusion/reaction time scales.We will also provide an alternative basis for the nullspace of the system with a barrier (Proposition 5.4).
To prove Theorem 2.1, we will apply concepts from linear perturbation theory and derive the SVD in the limit as diffusion becomes much faster than reactions.Specifically, we first consider the system at γ = 0, and derive a set of left singular vectors and singular values (i.e., the eigenvectors and eigenvalues of SS T when γ = 0).Because this system necessarily has repeating eigenvalues, the associated eigenvectors are not unique and are not necessarily continuous with respect to γ.However, we can apply results from Lemma 4.1 to find the unique orthonormal eigenprojection associated with each eigenvalue.Using these eigenprojections and perturbation theory results, we find the basis of eigenvectors that the system converges to continuously as γ → 0. For a review of the necessary concepts from perturbation theory that are used in the proof see Appendix A.
To prove Theorem 2.2, we show directly that the given equations are equivalent to the SVD.We also show that the SVD given by Theorem 2.2 is a simplified version of the SVD given by Theorem 2.1 (see Corollary 5.5).
5.1.The perturbed and unperturbed systems.The left singular vectors of S are given by the solutions to the following eigenvalue problem We will consider solutions to this eigenvalue problem in the limit as diffusion becomes much faster than reactions (i.e.γ → 0).To consider this in the context of perturbation theory, we rewrite the eigenvalue problem as follows (5.1) where now we are explicitly including the dependency of u i and λ i on γ.The unperturbed matrix is and the perturbation matrix is . Given appropriate choices for the eigenvectors u i (γ), the eigenvectors and eigenvalues will be continuous functions of γ in the neighborhood of γ = 0.
5.2.The eigenvalues and eigenprojections of the unperturbed system.In this section we provide an orthonormal eigendecomposition for the unperturbed matrix T (Lemma 5.1).We then use this eigendecomposition along with the results from Lemma 4.1 to find the unique orthonormal eigenprojections associated with each eigenvalue (Lemma 5.2).
Lemma 5.1.An orthonormal eigendecomposition of T is given as where The proof of this lemma is given in Supplemental Material D. From (5.10) it is immediately clear that T has repeating eigenvalues.This implies that the eigenvectors in the matrices given by (5.8) and (5.9) are not unique, and therefore, likely not the eigenvectors the system converges to as γ → 0.
Lemma 4.1 along with (5.10) imply that eigenvalues of T either repeat m + , m − , or m + m − times.Using this result and the set definitions defined in (3.7), we next identify each unique eigenvalue and find the associated orthonormal eigenprojection.
Lemma 5.2.The unique eigenvalues of T are contained in the following three sets (5.11) σ The corresponding unique orthonormal projections are, respectively, d ) T ⊗ D + j ∈ J (5.12) (5.13) where j 1 and j 2 are such that σ Proof.From Lemma 5.1, we see that every eigenvalue of T is contained in the sets defined by (5.11) and from Lemma 4.1 it follows that a given eigenvalue is only contained in one of the sets.Therefore, (5.11) contains the unique eigenvalues of T .
To find the orthonormal eigenprojection associated with each eigenvalue we will use the eigenvectors as defined in Lemma 5.1.Specifically, using the eigenvectors given by (5.8) and (5.9), we can use (A.4) in Appendix A to obtain the unique orthonormal eigenprojection.
For j ∈ J the m + eigenvectors associated with (σ (see Lemma 5.1).Therefore, the associated eigenprojection is (5.16) For j ∈ J C , there are m+m − eigenvectors associated with (σ . These eigenvectors are given by the columns of the following three matrices where j 1 and j 2 are as given by Lemma 4.1.Using the same logic as shown in (5.16), we have that the eigenprojection can be written as (5.17) Applying analogous logic for (σ (j) d1 ) 2 for j ∈ J 1 and (σ (j) d2 ) 2 for j ∈ J 2 leads to the eigenprojections P n1,j and P n2,j , respectively, as written in the claim.

The approximate left singular vectors of S.
We next use the eigenprojections given in Lemma 5.2 to derive the left singular vectors (i.e., eigenvectors of T ) that the system converges to continuously as γ → 0. This provides an approximate orthonormal basis for the column space and left null space of S as γ → 0.
Proposition 5.3.A complete set of left singular vectors of S and corresponding singular values is given by the columns/diagonal elements of the following matrices: .
Unlike Theorem 2.1, we are not identifying which singular vectors correspond to nonzero verse zero singular values.This is because the definitions provided in Proposition 5.3 allow for a direct comparison with the eigenprojections given by Lemma 5.2.However, Proposition 5.3 immediately gives the left singular vectors in Theorem 2.1.To see this note that where the matrices on the left of the equality represent those defined in Theorem 2.1 and the matrices on the right represent those defined in Proposition 5.3.The permutation matrices P 1 and P 2 are required to ensure the columns are in the correct order for comparison.Note that the singular values are related analogously.
Proof of Proposition 5.3.We will prove this result by considering the eigenprojections of T defined in Lemma 5.2.For each eigenprojection we calculate T (1) , given be (A.7), and its eigendecomposition.We will write (5.18) T (1) where • depends on the eigenvalue/eigenprojection we are considering and T (1) is given by (5.3).In the limit as γ → 0, the eigenvectors of (5.18) in the range of P• are equivalent to the left singular vectors.Additionally, the eigenvalues of (5.18) are used to find linear approximations of the singular values as shown by (A.8).First for j ∈ J , consider the eigenprojection given by (5.12).We have that (5.19) T (1) j = P j T (1) P j .
To find the eigendecomposition of (5.19), we apply Property C.1 given in Supplemental Material C to show that where recall S r+,j is given by (3.14).The eigenvectors of T (1) j in the range of P j are the columns of (5.20) Û1,j = u (j) and the corresponding eigenvalues are contained in (Σ

The eigenvectors of
T (1) n1,j and T (1) n2,j that are in the range of P n1,j and P n2,j are (5.24) and the corresponding eigenvalues are Σ 2 r1,− and Σ 2 r2,− , respectively.Recall the definition of S ri,− is given by (3.9).Again, using (A.8) this leads to the singular values given by Σ 2 and Σ 3 .
Finally, suppose j ∈ J C and consider the eigenprojection P j given by (5.13).Let j 1 and j 2 be as given by Lemma 4.1.For notational simplicity we will make the following substitutions We will also use the matrices B i for i = 1, ..., 5 given by (3.13).Using Property C.1 and C.2 given in Appendix C, the eigenvector relationship given by (4.1) in Lemma 4.1, and (5.18) we have that To obtain the eigendecomposition of T (1) j for j ∈ J C , we will suppose that the eigenvectors take the form and derive the values of v 1 , v 2 ∈ R m×1 .We have that (5.26) T (1) where we are using Property C.3 to calculate (uu T ⊗ B 1 )v.Therefore, for v to be an eigenvector of T (1) j the following smaller eigenvalue problem must hold where B is given by (3.12).This implies that [v 1 ; v 2 ] is equal to a left singular vector of S r , given by (3.10).Note that since [v 1 ; v 2 ] is a unit vector, v is also a unit vector and, thus, properly normalized.Only some of the singular vectors of S r result in eigenvectors v that are in the range of P j .Specifically, note the singular vectors contained in the columns of Ȗr,ex , see (3.15), result in eigenvectors that are not in the range of P j .To see this note that, when [v 1 , v 2 ] = Ȗr,ex , Coupled with the left singular vector relationship given by (4.1), this completes our derivation of the singular vectors contained in U 4 .Using (A.8), we obtain the singular values given by Σ 4 .

Right singular vectors.
Next we will approximate the right singular vectors of the system in the limit as γ → 0. To derive the right singular vectors that represent a basis for the row space, we use the following equation and the results from Proposition 5.3.For i = 1, ..., 5, (5.27) Vi = S T Ûi Σ−1 i Note that the equations for Ûi and Σi are given by Theorem 2.1, however their derivation is found in the proof to Proposition 5.3.Using this equation, we obtain the set of right singular vectors given by Theorem 2.1.
To complete the proof of Theorem 2.1, it remains to show that V defines a orthonormal basis for the nullspace that the system approaches as γ → 0. The complete proof of this is given in Supplemental Material D. Note that, an alternative asymptotic nullspace can be found.The nullspace given by Theorem 2.1 has the property that it is orthogonal for small values of γ and S V → 0 as γ → 0. It is possible to instead find a basis such that S V = 0 for small values of γ and the basis approaches orthogonal in the limit as γ → 0. The following lemma provides the equations for this alternative basis.
Proposition 5.4.The column vectors in the follow matrices span the nullspace of S as given by (2.1) and this basis is orthogonal in the limit as γ → 0: where where and Notice that only V1 and V3 have changed when compared to Theorem 2.1.
5.5.SVD for systems with spatially homogeneous reactions and diffusion.In the previous section we presented the approximate SVD for a system with a spatial barrier.Here, we will consider a specific scenario where there is no barrier and the reactions are the same across the domain.In terms of the previous notation, this is equivalent to setting m + = m, m − = 0 and S r1 = S r2 .Under these conditions, we will show that the SVD reduces to a simplified form (Corollary 5.5 and 5.6) and becomes exact for all values of γ, i.e., prove Theorem 2.2.Below we set S r = S r1 and refer to the singular value decomposition of S r using the notation given in (2.2).
First note, that under these conditions the SVDs of the stoichiometry-like matrices are simplified.We have that the SVD of S r+,j is (5.30) This result is shown by considering the equations for S r and S r+,j as given by (3.10) and (3.14), respectively.Using (5.30) and (5.31), we next show that the SVD given by Theorem 2.1 reduces to a simplified form.
Corollary 5.5.The left singular vectors given by Theorem 2.1 reduce to the columns of following matrix and the singular values reduce to the diagonal of Proof.To prove this corollary we will examine the SVD for a system with a barrier.Specifically, we consider the left singular vectors and values as written in Proposition 5.3.We will show that these vectors and values reduce to the singular vectors and values given in the corollary statement First, note that U 2 and U 3 are empty matrices since M − is an empty set.For U 1,j , j ∈ J , using that I M+ m = I m .and the SVD given by 5.30 we have that, (5.34) Û1,j = u For U 4 using the SVD given by 5.31 and that q r = q r and qr = m − q r , we have that .
Putting these results together we obtain the set of left singular vectors and singular values given by the corollary statement.
Corollary 5.6.The right singular vectors of S given by Theorem 2.1 reduce to the following for the simplified system.
Proof.For the right singular vectors we will use the equations as given in Theorem 2.1.We will show the proof for V and note that the proof for V follows analogously.Note that V2 , V3 , and V4 are empty.
For V1,j , using the SVD given by (5.30), we have that For V5 and V6 , using the SVD given by (5.31), we have that Putting this together and rearranging columns we obtain the equation for V given in the corollary statement.
Finally, we will prove the main result that the SVD of the simplified S is valid for all values of γ.
Proof of Theorem 2.2.To show that (2.9) is the SVD of S as given by (2.8), it is suffices to show that, U and V are orthogonal matrices and S = Û Σ V T .
We first show that U and V are orthogonal matrices.Recall that U ∈ R nm×nm where U = U d ⊗ U r .We have that It follows that U is orthogonal.
For V , we have that where recall that Σ2 = Σ2 d ⊕ γ Σr

2
. It can similarly be shown that V T V = 0, V V T = 0, and V T V = I.
Next, we will show that S = Û Σ V T .
Therefore, Theorem 2.1 gives the SVD of S at all values of γ.
5.6.Error analysis for example system.In this section we present an error analysis for the approximate SVD of an example stoichiometry matrix.We demonstrate numerically that the approximate SVD presented in Theorem 2.1 converges to the true SVD in the limit as γ → 0. We will consider a simplified set of equations that describes part of the Calvin Cycle in cyanobacteria.Specifically, cyanobacteria have cellular compartments called carboxysomes that serve to concentrate carbon within the cell [12].It is thought that this compartmentalization increases the amount of carbon fixation and decreases the flux through the competing photorespiration reaction.This is an example of the type of system that could, in the future, be investigated more thoroughly with the approach presented here.We consider a system with n = 8 compartments, where n 1 = 2 and n 2 = 6.The first subregion in the domain represents the carboxysome and the second region represents the cytoplasm.We will consider the scenario of Mixed boundary conditions where fluxes are allowed only into the right side of the domain (i.e., into the cytoplasm region).Biologically, this scenario could represent a radially symmetric region in the cell centered on a carboxysome.The species in this system, as ordered in the stoichiometry matrix, are bicarbonate (HCO 3 -), Ribulose 1,5-bisphosphate (RuBP), carbon dioxide (CO 2 ), 3-phosphogylcerate (3 PGA), Oxygen (O 2 ), and 2phosphoglycolate (2 PG).The reactions are given as (5.41) It is known that O 2 and CO 2 cannot diffuse into the carboxysome [7,9].Therefore we set M − = {3, 5} and M + = {1, 2, 4, 6}.
Given that O 2 is not present in the carboxysome, we know that only R1 and R3 occur in the first subregion.This leads to the following reaction-only stoichiometry matrices in the first and second region, respectively, (5.42) Using the defined parameters, we applied the equations in Theorem 2.1 at multiple values of γ and compared the results to the numerical SVD in MATLAB (Figure 1).As expected we find that the error approaches zero as γ → 0. In this comparison the singular vectors/values are sorted by the magnitude of the singular value.Singular values between the numerical and approximate SVD (and hence singular vectors) are paired by finding those that are closest to each other in size.Note that in the error analysis in Figure 1, we only consider the nonzero singular values and corresponding singular vectors.Similar results are observed for the right and left null space (e.g., S V → 0 as γ → 0).

Discussion
. In this paper we derived the approximate SVD of the stoichiometry matrix for a one dimensional discrete reaction-diffusion system partitioned into two subregions.Between these two subregions only certain species are allowed to diffuse.We additionally presented the exact SVD in the scenario where diffusion is allowed freely throughout the domain.This work provides a framework that can be applied and expanded upon to examine a variety of reaction-diffusion scenarios.For example, we hypothesize that in more complex scenarios (e.g., species-dependent boundary conditions) a Kronecker product formulation can still be used to write the SVD.Additionally, the formulas given by Theorem 2.1 and 2.2 allow for future analysis looking at the effects of spatial properties in compartmentalized systems.
Computationally, the results of Theorem 2.1 and 2.2 allow for the efficient estimation of the SVD of the RD stoichiometry matrix.Importantly, the approximate SVD is fully determined by the SVDs of smaller matrices.Either the SVD of these smaller matrices is known analytically or the dimension of the matrices is independent of the number of spatial compartments in the system.For example, consider a system with m species, p reactions, and either n or 2n compartments.The diffusion-only stoichiometry matrices for this system are known analytically.The other matrices that the SVD depends on have dimensions proportional to m and/or p. Notably the total number of required smaller matrix decompositions will increase linearly with the number of compartments.6.1.Intuition for SVD results.The approximate SVD for a system with a barrier provides intuition for how the system's structure influences dynamical and steady state properties.As written in Theorem 2.1, we have partitioned the singular vectors into multiple sets, which we will refer to as eigenreaction sets.For example, the singular vectors in Û1,j and V1,j and the singular values in Σ1,j for j ∈ J represent the first eigenreaction set.Recall that the SVD defines eigenreactions, which represent decoupled linear combinations of species that are moved by a linear combinations of fluxes, e.g., (6.1) where w is the species concentration vector and f is the vector of fluxes, e.g., see (1.1).Each eigenreaction set describes the movement of species with similar diffusive properties.The first eigenreaction set describes the movement of species that are able to diffuse across the barrier.The second, third and fourth eigenreaction sets all describe the movement of species that are unable to diffuse across the barrier.The second eigenreaction set includes both reactive and diffusive movement in the first subregion, the third eigenreaction set describes only reactive movement in the first subregion, and, finally, the fourth eigenreaction set describes reactive and diffusive movement in the second subregion.Recall that the third eigenreaction set is only nonempty for Mixed boundary conditions.The fifth and sixth eigenreaction sets are unique in that they describe the movement of all species in the system.This movement is coupled due to the repeating singular values in the diffusion-only stoichiometry matrices (i.e., S d , S d1 , and S d2 ).This demonstrates that even in the regime where diffusion is much faster than reactions, there is still a coupling between species with different diffusive processes.
The basis for the nullspace of S is also partition into multiple sets, given by Theorem 2.1 (i.e., Vi for i = 1, ..., 5).We will refer to these as steady-state flux sets, since they represent fluxes that can exist under steady-state conditions.The first and third steady-state flux sets include only reactive fluxes in the first subregion and throughout the domain, respectively.Note that the first steady-state flux set is only nonempty for Mixed boundary conditions, whereas the third steady-state flux set is only nonempty for Zero Flux boundary conditions.The second and fourth steadystate flux sets represent reactive and diffusive flux combinations in the first subregion and second subregion, respectively.Finally the fifth steady-state flux set represents fluxes that are in the nullspace of S due to their infeasibility.That is, for a given dynamical system, these fluxes will never contribute since they define fluxes across barriers/boundaries that are not allowed.

Conclusion and Future Work.
To find the SVD of the RD stoichiometry matrix, we first used linear perturbation theory to calculate the left singular vectors and values.We then used the resulting vectors and values to find the right singular vectors.An alternative approach would be to instead derive the right singular vectors directly using perturbation theory.Although this approach may provide additional insight into the system properties, it is slightly more complex as it involves additional terms in the expansions used in the perturbation analysis.Therefore, this analysis is the topic of future research.
The key assumption used to derive the approximate SVD is that diffusion is much faster than the reactions.Whether this is a valid assumption depends on the specific biological system under consideration.Indeed, the relative time-scales of diffusion and reactions in biological systems can very greatly and is a complex topic [13].A similar approach, as presented in this paper, could be applied to derive the approximate SVD in a system where reactions occur much faster than diffusion.That is, we would instead consider the perturbation problem in the limit as the diffusive term of (5.1) goes to zero.Rigorously showing whether the approach applied here could work in this alternative case is a topic of future research.
Our motivation in deriving the approximate SVD in terms of reduced systems is to gain insight into how including spatial barriers and diffusion impact a biological system.The SVD for the reaction-only system has provided valuable insight in comparing genome-scale metabolic networks [4] and finding connections between biochemical processes [10].By including spatial parameters, the work presented here provides tools that computational biologists can apply to understand how reactive processes are coupled across space.
The SVD will depend on the following constants for j = 1, ..., n, Next we define the left singular vectors that correspond to the column space and left nullspace.For the left singular vectors in the column space, the ith element of the jth left singular vector is, for i = 1, ..., n and j = 1, ..., q d , (B.4) u The left nullspace is only nonempty for Zero flux boundary conditions and we have that For the right singular vectors, the ith element of the jth right singular vector associated with nonzero singular values is, for i = 1, ..., n + 1 and j = 1, ..., q d , (B. 6) Open.
For the right singular vectors in the nullspace of S d , for Zero Flux boundary conditions, we have that = e n+1 (B.7) where e i represents the vector with zeros and a one at the ith index.For Mixed, Mixed-Alt, and Open boundary conditions, we have Mixed-Alt where 1 is a vector of ones.
Finally, the jth singular value, for each of the boundary conditions, is 2 sin (a n,j ) Zero Flux 2 sin (b n,j ) Mixed and Mixed-Alt 2 sin (c n,j ) Open.
Appendix C. Kronecker product formulas.
In this section we provide some Kronecker product relations that are needed to prove Theorem 2.1.We omit the proof of these properties but note that they can be shown through a series of algebraic manipulations.
Let A 1 , A 2 , B 1 , B 2 be square matrices of the same size.Then (C.2) Property C.2. Suppose u 1 and u 2 are unit vectors and a i ∈ R for i = 1, .., 8. Let A 1 , A 2 , B 1 , B 2 be square matrices of the same size.Then Property C.3.Suppose that u 1 and u 2 are unit column vectors and a i ∈ R for i = 1, .., 6 such that and a 4 = a 2 a 6 a 5 .
Let A 1 ∈ R m×n and v 1 , v 2 ∈ R m×1 .We then have that Appendix D. Proofs.This section contains supplemental proofs for the results presented in Section 4 and 5.We first provide the proof of Lemma 4.1, which provides a relationship for the eigenvalues and eigenvectors of S d , S d1 and S d2 .
Proof of Lemma 4.1.We will consider the three possible boundary conditions independently.
which completes the proof.
Next, we prove Lemma 5.1, which provides an eigendecomposition of SS T when γ = 0, recall S is given by (2.1).This is equivalent to the non-unique eigendecomposition of the unperturbed matrix T , given by (5.2).
Proof of Lemma 5.1.First note that T ∈ R nm×nm and, as needed, the number of eigenvectors defined is nm + + n 1 m − + n 2 m − = nm.
We will show that the matrices QT,1 , QT,2 , and QT,3 contain eigenvectors of T and ΣQ T ,1 , ΣQ T ,2 , and ΣQ T ,3 contain the corresponding nonzero eigenvalues.First, considering QT,1 , we have that We leave it as an exercise to show that the eigenvectors and nullspace basis vectors are orthogonal.
Next, we will show that Theorem 2.1 provides an approximate basis for the nullspace of S (i.e., V ), where the basis is orthogonal at small gamma and satisfies S V = 0 in the limit as γ goes to zero.We will also provide the proof to Proposition 5.4, which gives an exact basis for the nullspace of S that is orthogonal in the limit as γ → 0.
Proof of Theorem 2.1 (nullspace).The dimension of the nullspace of S is given by Lemma 4.2.Notice that this dimension matches the number of columns in V as defined in Theorem 2.1.Specifically, for the five matrices that compose V , i.e.Vi for i = 1, ..., 5 the number of columns is and by inspection we see that the number of columns is equivlanet to the value of q given by Lemma 4.2.
We leave it as an exercise to show that all the vectors defined in these matrices are orthonormal.
To show that the vectors are in the nullspace of S, write S as follows where S d,s1 represents the first n 1 rows of S d and S d,s2 represents the last n 2 rows of S d .Suppose a vector in the nullspace can be written as where Σ is a diagonal matrix.Multiplying S by v, we obtain the following two equations that must be satisfied It is straightforward to show that the vectors given by the claim satisfy these equations in the limit as γ → 0. In fact, for V2 , V4 and V6 the equations are satisfied at small gamma.Below we will show the logic for V2 .We leave it as an exercise to verify these results for V4 and V6 .Additionally, it is trivial to show that as γ → 0, V1 and V3 satisfy the condtions since, in this case, v 5 , v 6 = 0.
For V2 we have that Proof of Proposition 5.4.The proof to this proposition closely follows the proof given for the nullspace presented in Theorem 2.1.In addition to the logic of this proof we need to show that the basis vectors that differ (i.e., those in V1 and V3 ) satisfy the two conditions given in D.7 and D.8 at small values of γ.We will show the logic for V3 and leave it as an exercise to show that V1 satisfies the conditions.
For V3 , when considering the conditions given by D.7 and D.8, we have that Ȗd,s2 Ûd,s1 Û T d,s2 and C 2 /C 1 = n 2 /n 1 .We leave it as an exercise to show that all the vectors defined in the columns of V are linearly independent and that, in the limit as γ → 0, they become orthogonal.

Specifically, b 1
= 0, b 2 = 0 implies zero flux boundary conditions, b 1 = 0, b 2 = 1 implies Mixed boundary conditions and b 1 = 1, b 2 = 1 implies Open boundary conditions.The diffusion only-stoichiometry matrices S d1 and S d2 are defined similarly.However, for S d1 the value of b 2 is replaced by zero and for S d2 the value of b 1 is replaced by zero.The n rows of S d corresponds to the species in each of the n compartments, and the n + 1 columns corresponds to the flux across each of the n − 1 interior edges as well as the two boundaries at either end of the domain.Using the following matrix, S r1,+ := (S r1 ) M+ S r1,− := (S r1 ) M− S r2,+ := (S r2 ) M+ S r2,− := (S r2 ) M− where (A) B represents the rows of A that are in the index-set B.

T
where ΛT contains the nonzero eigenvalues of T and QT = QT,1 QT,2 QT,

Fig. 1 .
Fig. 1.Error of approximate SVD at decreasing values of γ.The 'num' subscript refers to singular vectors/values obtained numerically using MATLAB.

Next for QT, 2 ,
we have thatT QT,2 = S d S T d ⊗ D + + S d S T d − HH T ⊗ D − Ûd1 0 ⊗ I M− m = (S d S T d − HH T ) Ûd1 0 ⊗ I M− m = Ûd1 0 Σ2 d1 ⊗ I M− m = QT,2 Σ2 d1 ⊗ I m− .Finally, for QT,3 we have thatT QT,3 = S d S T d ⊗ D + + S d S T d − HH T ⊗ D −It can analogously be shown that QT,1 , QT,2 , and QT,3 represent the nullspace of T .

1
The Kronecker product of A ∈ R m×n and B ∈ R p×r is a mp + nr block matrix where

Table 1
Description of constants, index sets, and stoichiometry/stoichiometry-like matrices used to define the SVD of the RD stoichiometry matrix.