# Chemical Reaction Network Theory

In this post, we will introduce one of the applications of Numerical Algebraic Geometry, which is called the Chemical Reaction Network Theory. Many objects in Chemistry and Biology, for example cells, molecules and so on, make decisions by themselves. Because of these decisions, they decided to emerge each other, move around, make new material and even die! In other words, many causes in natural environment make them react, and that reaction is the decision that they make. Now, let’s elaborate on how Numerical Algebraic Geometry is applied on analysing chemical reactions.

In order to explain readily, we will use examples. Consider the following two chemical reactions.

In here, we call $k_1,\dots, k_6$ as rate constants which determines the frequency of the reaction. With this family of chemical reactions, we consider the triple $\mathcal{N}$ which is called the network. The network $\mathcal{N}$ is determined by $3$ components, $\mathcal{S}$ specie, $\mathcal{C}$ complex, and $\mathcal{R}$ reaction (hence we can write $\mathcal{N}=(\mathcal{S},\mathcal{C},\mathcal{R})$ as a network. Specie is defined by the set of components which are involved in the reactions. Therefore, we have $\mathcal{S}=\{A,B,C,D,E\}$ in the example above. Complex is the set of vertices in the chemical reaction diagram, which is $\mathcal{C}=\{A, 2B, A+C, D, B+E\}$. Finally, reaction is the set of edges.

From now on, we will observe the Mass-Action Kinetics only. The Mass-Action Kinetic is the family of chemical reactions which the reaction rate is proportional to the product of the concentration of species.

Then, according to the definition of the Mass-Action Kinetics, we can construct the ODE system of the concentration of species. First, let $x_A, x_B, x_C, x_D, x_E$ be the concentration of $A,B,C,D,E$ respectively. Then, have the following ODEs for the concentration of species:

$\dot{x_A}=k_1x_B^2-k_2x_A+k_3x_D-k_4x_Ax_C+k_5x_Bx_E\\ \dot{x_B}=-2k_1x_B^2+2k_1x_A-k_5x_Bx_E+k_6x_D\\\vdots$

If you are too tired to write down, then you can choose the following way. Let $y\in \mathbb{Z}^{n\times d}$ where $|\mathcal{S}|=d$ and $|\mathcal{C}|=n$ whose $i$-th row is the exponent vector of the monomial corresponding to the $i$-th complex. For example, for $\mathcal{C}=\{A,2B,A+C,D,B+E\}$, we have $x_A,x_B^2,x_Ax_C,x_D,x_Bx_E$ and so, $y=\left[\begin{array}{ccccc} 1&0&0&0&0\\0&2&0&0&0\\1&0&1&0&0\\0&0&0&1&0\\0&1&0&0&1 \end{array}\right]$. In addition, let $\Psi(x)$ be the vector of the monomials. Thus, we have $\Psi(x)=(x_A,x_B^2,x_Ax_C,x_D,x_Bx_E)$. Finally, we define the matrix $A$ as the negative laplacian matrix of the graph. For example, we have $A=\left[\begin{array}{ccccc} -k_2&k_2&0&0&0\\k_1&-k_1&0&0&0\\0&0&-k_4&k_4&0\\0&0&k_3&-k_3-k_6&k_6\\0&0&k_5&0&-k_5 \end{array}\right]$ in our example. Then, we have the same system of ODEs $\frac{dx}{dt}=\Psi\cdot A \cdot y$. By letting $\frac{dx}{dt}=0$, we can apply the concept of the variety on polynomial equations. We call that variety as a steady state variety.

Moreover, by using this network, we can define the reaction vector. which is the difference of the coefficient vectors that involved in each reaction. For example, if we have the reaction $A\rightarrow B$, then we have the corresponding reaction vector $\left[\begin{array}{c} -1\\2\\0\\0\\0 \end{array}\right]=\left[\begin{array}{c} 0\\2\\0\\0\\0 \end{array}\right]-\left[\begin{array}{c} 1\\0\\0\\0\\0 \end{array}\right]$ (consider $2B-A$). The span of reaction vectors $S$ is called the stoichiometric subspace. If we have the initial condition $c_0$, then we have the $S+c_0$ which have the same dimension and degree with $S$. The stoichiometric subspace have the meaning that whatever the time is the amount of complexes are contained in $S$ because every element in $S$ is spanned by the results of reactions.

With these knowledge, we may have several questions to ask. Making computation code (especially on Macaulay2) to compute the ODE system by using network, simplifying system of ODEs by deleting one of reactions (that is, deleting the variable $k_i$). If we simplify the network, then how degree and dimension are changed? If we delete $k_i$ by converging $k_i$ to $0$ with the method of homotopy continuation, then what is the meaning of quadratic turning point in that case?