# Mass Action Law

## 1. Introduction

Let's assume a system of interacting species (chemicals, cells, genes, proteins, organisms...). The ```Law of Mass Action``` states that the change in the amount of a given species is proportional to the product of the amount of the interacting substances.  Lets define some species `A`, `B`, `C` and `D` that interact between them, so the amount of each one is allowed to change over time until the reaction comes to an end. 

Some possible schemes of interaction are, for instance:
$$
\begin{align*}
A + B  &\overset{k_1}{\longrightarrow} C  \tag{1} \\
2A + B  &\overset{k_2}{\longrightarrow} C  \tag{2} \\ 
C  &\overset{k_3}{\longrightarrow} A + B   \tag{3} \\ 
AB + C  &\overset{k_4}{\longrightarrow} AC + B   \tag{4} \\ 
AB + CD  &\overset{k_5}{\longrightarrow} AC + BD   \tag{5} \\ 
\end{align*}$$

The speed of a given reaction (i.e., the change in the amount of a given reactant) is proportional to the amount of collisions (the amount of times that two given reactants crash into each other in the solution), and this rate of collisions is proportional to the amount of reactants in the solution. For instance, for reaction 1


```julia
using Plots 
p1 =scatter(fill(rand(10),1), fill(rand(10),1), fill(rand(10),1), m=(10, 0.8, :blues),framestyle=:box, title="few collisions",label="A")
p1 =scatter!(fill(rand(10),1), fill(rand(10),1), fill(rand(10),1), m=(10, 0.8, :reds),framestyle=:box, title="few collisions",label="B")
p2 =scatter(fill(rand(50),1), fill(rand(50),1), fill(rand(50),1), m=(10, 0.8, :blues),framestyle=:box, title="Many collisions",label="A")
p2 =scatter!(fill(rand(50),1), fill(rand(50),1), fill(rand(50),1), m=(10, 0.8, :reds),framestyle=:box, title="Many collisions",label="B")
plot(p1, p2)
```




    

    



Therefore, the speed of each reaction above can be written as:

$$
\begin{align*}
speed~of~reaction~1~&\propto A \cdot B \tag{6} \\ 
speed~of~reaction~2~ &\propto A \cdot A \cdot B \tag{7} \\ 
speed~of~reaction~3 ~ &\propto C \tag{7} \\ 
speed ~of ~reaction~4 ~&\propto AB \cdot B    \tag{8} \\ 
speed ~of ~reaction~5 ~&\propto AB \cdot CD  \tag{9} \\ 
\end{align*}$$

The proportionality constant is the value of the reaction rates $k_i$, which is a measure of how many of the collisions are effecive, so:

$$
\begin{align*}
speed~of~reaction~1~&= k_1 \cdot A \cdot B \tag{10} \\ 
speed~of~reaction~2~ &= k_2 \cdot A \cdot A \cdot B \tag{11} \\ 
speed~of~reaction~3 ~ &= k_3 \cdot C \tag{12} \\ 
speed ~of ~reaction~4 ~&= k_4 \cdot AB \cdot B    \tag{13} \\ 
speed ~of ~reaction~5 ~&= k_5 \cdot AB \cdot CD  \tag{14} 
\end{align*}$$

The equilibrium is reached when the speed is zero, i.e., when the reaction is finished because one of the reactants has been fully consumed. This statement is correct in irreversible reactions. On the contrary, in reversible reactions, the equilibrium is a more interesting concept. For instance, let's asumme a very simple equilibrium between two forms of a chemical:

$$
\begin{align*}
2 NO_2   &\overset{k_1}{\underset{k_2}{\longleftrightarrow}} N_{2}O_4 \tag{15} \\
\end{align*}
$$

This system is actually composed of two different reactions, 

$$
\begin{align*}
2 NO_2   &\overset{k_1}{\longrightarrow}  N_{2}O_4 \tag{16} \\
N_{2}O_4  &\overset{k_2}{\longrightarrow} 2 NO_2 \tag{17}
\end{align*}
$$

and the speed of each reactant is

$$
\begin{align*}
speed~of~reaction~16~&= k_1 \cdot NO_2 \cdot NO_2 = k_1 \cdot NO_{2}^{2} \tag{18} \\ 
speed~of~reaction~17~ &= k_2 \cdot N_{2}O_4 \tag{19}
\end{align*}$$


At room temperature, it is impossible to have pure $NO_2$ or  $N_{2}O_4$, so the equilibrium cannot be defined as the point when no reactions are taking place, or the point where one of the rectats has been fully consumed. In these conditions, the equilibrium is better defined as the point in time when there is no change in the amount of the reactans. In other words, the equilibirum occurs when the speed of the two reactions is equal. 
$$
\begin{align*}
k_1 \cdot NO_{2}^{2} = k_2 \cdot N_{2}O_4  \tag{20} 
\end{align*}$$

Thererefore, in conditions of equilibrum in reversible reactions, the ratio $\frac{[NO_2]^2}{[N_{2}O_4]}$ is constant and proportional to the ratios between the reaction rates. Interestingly, this ratio is allways the same independently on the amount of reactats that you added in initially in the system. This ratio depends on the temperature, since the reactions rates depend on the temperature (the number of collisions depens on the temperature). Such a law of nature is called the Law of Mass Action or mass action law. Although the Law of Mass Action was initially derived in the context of chemical reactions, it can be used also to describe the equilibirum of any system of interacting entities. We can define an equilibrium constant $K_{eq}$ as the ratio between the final concentrations of the two reactants $\frac{[NO_2]^2}{[N_{2}O_4]}$. The value of the equilibrium constant is then 
$$
\begin{align*}
K_{eq}=\frac{k_2}{k_1} \tag{21} 
\end{align*}$$


```julia
using Plots 
NO2= 0.1*LinRange(0,1,100)
K_eq=0.3
plot(NO2,K_eq*NO2.^2)
title!("Equilibrium concentrations")
xlabel!("NO_2 [Mol]")
ylabel!("N_2O_4 [Mol]")
```




    

    



## 2.General Formulation of the Mass Action Law
We introduce the mass action law by using a general chemical reaction equation in which species `A` and `B` react to give species `C` and `D`

$$
\begin{align*}
aA + bB  &\overset{k_1}{\underset{k_2}{\longleftrightarrow}} cC + dD \tag{22} \\ 
\end{align*}$$

where `a`, `b`, `c`, `d` are the stoichiometric coefficients for a balanced interaction. At any instant in time during the course of the reaction, we can define a relation between the values of each concentration:

$$
Q= \frac{[C]^c[D]^d}{[A]^a[B]^b} \tag{23}
$$
where `Q` is defined as the reaction quotient, and measures the relative amounts of the interacting species present during a reaction at a particular point in time. As the reaction takes place and the system moves towards equilibrium, the value of `Q` gradually approaches to the equilibrium constant $K_{eq}$. The Mass Action Law states that if the system is at its equilibrium, then the following ratio is a constant:

$$
K_{eq}= \frac{[C]^c[D]^d}{[A]^a[B]^b} \tag{24}
$$

The values of `[A]`, `[B]`, `[C]` and `[D]` represent the amount of each species at equilibirum. If this amount is given as a concentration, lets say in moles/volume, `M`, the units of the equilibrium constant $K_{eq}$ are:

$$
[K_{eq}]=M^{c+d-a-b} \tag{25}
$$

For instance, lets consider the following reaction

$$NaCO_3 + CaCl_2 \overset{k_1}{\underset{k_2}{\longleftrightarrow}} CaCO_3 + 2 \cdot NaCl \tag{26}\\  $$
The equilibrium constant is

$$
K_{eq}= \frac{[CaCO_3][NaCl]^2}{[NaCO_3][CaCl_2]} \tag{27}\\
$$

and the units are

$$
[K_{eq}]=M^{1+2-1-1}=M  \tag{28}\\
$$

If we can measure the final concentration at equilibrium of each species, for instance:


```julia
NaCO_3=2;
CaCl_2=0.5;
CaCO_3=2;
NaCl=1.2;
```

Then the equilibrium constant can be calculated as:


```julia
K_eq=(CaCO_3*NaCl^2)/(NaCO_3*CaCl_2)
```




    2.88



In conclusion, the ``Law of Mass Action`` says that the ratio between the concentrations of the reactants at equilibrium depends only on the ratio between the kinetic constants, and it is independent on the initial concentration of the reactants. To illustrate that, we will solve the reaction and plot the value of ``Q(t)`` 


```julia
using DifferentialEquations
using ParameterizedFunctions
```


```julia
simpleODE3! = @ode_def abetterway begin
  da = -k1 * a * b + k2 * c * d^2
  db = -k1 * a * b + k2 * c * d^2
  dc = k1 * a * b - k2 * c * d^2
  dd = 2 * k1 * a * b - 2 * k2 * c * d^2
    end k1 k2
```




    (::abetterway{var"#19#23",var"#20#24",var"#21#25",Nothing,Nothing,var"#22#26",Expr,Expr}) (generic function with 2 methods)




```julia
k1=2.3e0;  # units 1/(Ms)
k2=2.5e0;  # units 1/(M M s)

p = (k1,k2)

a₀=0.02; # units (M)
b₀=0.01; # units (M)
c₀=0; # units (M)
d₀=0; # units (M)

u₀=[a₀,b₀,c₀,d₀];
```


```julia
tspan = (0.0,300.0)
prob = ODEProblem(simpleODE3!,u₀,tspan,p)
sol = solve(prob)
plot(sol,title= "Solution of differential equations",xlabel="Time [s]",ylabel="Concentration [M]")
```




    

    




```julia
plot(sol.t,(sol[3,:].*sol[4,:].^2)./(sol[1,:].*sol[2,:]),label="Q",title= "Plot of the quotient coefficient",xlabel="Time [s]",ylabel="Q [M]")
```




    

    



Now we solve the system again but with different initial conditions. We see that, although the solution is very different, the value of $K_eq$ (the value of $Q$ at equilibrium) is maintained. 


```julia
k1=2.3e0;  # units 1/(Ms)
k2=2.5e0;  # units 1/(M M s)

p = (k1,k2)

a₀=0.05; # units (M)
b₀=0.01; # units (M)
c₀=0; # units (M)
d₀=0; # units (M)

u₀=[a₀,b₀,c₀,d₀];

tspan = (0.0,300.0)
prob = ODEProblem(simpleODE3!,u₀,tspan,p)
sol = solve(prob)
plot(sol,title= "Solution of the differential equations",xlabel="Time [s]",ylabel="Concentration [M]")
```




    

    




```julia
plot(sol.t,(sol[3,:].*sol[4,:].^2)./(sol[1,:].*sol[2,:]),label="Q",title= "Plot of the quotient coefficient",xlabel="Time [s]",ylabel="Q [M]")
```




    

    



## 3. Matrix Formulation for systems of interactions

Sometimes, it is convenient to write down the Mass Action Law as a vector-matrix formulation. This will be useful when dealing with large numbers of species and or reactions. Consider the familiar reaction:

$$ 2H_2 + O_2 \overset{k}{\longrightarrow}  2 H_2O \tag{29}$$

Eq. 12 can be rewritten in a more general form as
 
 $$ 2 X_1 + X_2 \overset{k}{\longrightarrow}  2 X_3 \tag{30}$$

where $X_1$, $X_2$, and $X_3$ are the species $H_2$, $O_2$, and $H_2O$,
 
The first step to write the system as a matrix form, we place all reactants at both sides of the reaction.

  $$ 2 X_1 + X_2  + 0 X_3 \overset{k}{\longrightarrow}  0 X_1 +  0 X_2  + 2 X_3 \tag{31}$$
 
The values in front of each number are the stoichimetric coefficients, which state the number of molecules that participate in each reaction. We can write these stoichiometric coefficiens as
 
$$ A_1 X_1 + A_2 X_2  + A_3 X_3 \overset{k}{\longrightarrow}  B_1 X_1 +  B_2 X_2  + B_3 X_3 \tag{32}$$

$A_1$ = 2, $A_2$ = 1, $A_3$ = 0, $B_1$ = 0, $B_2$ = 0, and $B_3$ = 2.

In compact form, the reaction can be written as:

$$
\sum_{j=1}^{3} A_j X_j \overset{k}{\longrightarrow}  \sum_{j=1}^{3} B_j X_j \tag{33}
$$

or in vector form:$ A X \overset{k}{\longrightarrow} B X  \tag{34}$

where `X` is the state vector (i.e, the vector with all reactants), and `A` and `B` are the matrices of stoichiometric coefficients 

$$\begin{align*}
X&=[X_1,X_2,X_3]^T=[H_2,O_2,H_2O]^T \tag{35}\\
A&=[A_1,A_2,A_3]=[2,1,0] \tag{36}\\
B&=[B_1,B_2,B_3]=[0,0,2] \tag{37}\\
\end{align*}$$

To learn how to write reversible reactions as a matrix form, let's practice with the following reaction:

$$NaCO_3 + CaCl_2 \overset{k_1}{\underset{k_2}{\longleftrightarrow}} CaCO_3 + 2 \cdot NaCl \tag{38}\\  $$

We will decouple the reversible reaction into two irreversible steps:

$$\begin{align*}
NaCO_3 + CaCl_2 \overset{k_1}{\longrightarrow} CaCO_3 + 2 \cdot NaCl \tag{39}\\ 
CaCO_3 + 2 \cdot NaCl \overset{k_2}{\longrightarrow} NaCO_3 + CaCl_2 \tag{40}
\end{align*}$$

Using the same matrix notation as in Eq. 13: 

$$\begin{align*}
X_1 + X_2 \overset{k_1}{\longrightarrow} X_3 + 2 \cdot X_4 \tag{41}\\ 
X_3 + 2 \cdot X_4 \overset{k_2}{\longrightarrow} X_1 + X_2 \tag{42}
\end{align*}$$


where $X_1 = Na_2CO_3$, $X_2  = Ca_2Cl$,$X_3  = CaCo_3$, $X_4 = NaCl$. Now we put all reactants at both sides of the reaction.  

$$\begin{align*}
1 \cdot X_1 + 1 \cdot X_2  + 0 \cdot X_3  + 0 \cdot X_4 & \overset{k_1}  \longrightarrow  0 \cdot X_1+ 0 \cdot X_2 + 1 \cdot X_3 + 2\cdot  X_4 \tag{43}\\
 0 \cdot X_1+ 0 \cdot X_2 +1 \cdot X_3 + 2\cdot X_4 &\overset{k_2} \longrightarrow 1 \cdot X_1 + 1\cdot  X_2 + 0 \cdot X_3 + 0 \cdot X_4 \tag{44}
\end{align*}$$

in a more general form:
$$\begin{align*}
A_{11} X_1 + A_{12} X_2  + A_{13} + A_{14} X_4 &\overset{k_1}  \longrightarrow  B_{11}X_1+  B_{12} X_2 +  B_{13} X_3 + B_{14} X_4 \tag{45}\\
A_{21} X_1+ A_{22} X_2 +A_{23} X_3 + A_{24} X_4 & \overset{k_2} \longrightarrow B_{21}  X_1 + B_{22}  X_2 + B_{23} X_3 + B_{24} X_4 \tag{46}
\end{align*}$$

The numbers $A_{ij}$ are the stoichiometric coefficientes of the reactants, $B_{ij}$ are the stoichiometric coefficientes of the products.  Using the notation of Eq 16, these values are now the components of matrices `A` and `B`:

$$
A=\begin{bmatrix}
 1& 1 & 0 & 0  \\ 
 0 & 0 &1   &2  
\end{bmatrix} ;
B=\begin{bmatrix}
0 & 0 &1   &2 \tag{47}\\ 
1 & 1 &0   &0  
\end{bmatrix}
\\
$$

The numbers of rows $i= 1...r$, corresponds to the numbert of reactions `r` involved. The number of columns $j= 1...s$ corresponds to the number of species involved `s`. In conclusion, when we want to write the mass action law for `r` reactions, `k` is now also a vector with components:

$$k=[k_1,k_2,...k_r]^T \tag{48}$$


altogether we have 


$$
\sum_{j=1}^{s} A_{ij} X_j \overset{k}{\longrightarrow}  \sum_{j=1}^{s} B_{ij} X_j,\,\,\,\, \,\,\,\,\,\,\,\, \,\,\,\,  i=1...r\tag{49}
$$



### Exercise: 
White a computer program that consists on a function that takes as input the matrices A and B and returns the value $(B-A)'$

## 4. Order of reaction and units of k

The Order of a Reaction refers to the power dependence of the rate on the concentration of each reactant. In brief, the order of the reaction indicates the correlation of its velocity to the amount of reactants.
- For a ``zero-order reaction``, the rate does not depend on the concentration of any species.
- For a ``first-order reaction``, the rate is dependent on the concentration of a single species.
- For a ``second-order reaction``, the rate is dependent on the square of the concentration of a single reactant, or two reactants.

In Mass-Action kinetics, we can define the reaction order of a given reaction $i =1…r$ as:

$$
\sum_{j=1}^{s} A_{ij} \tag{50}
$$

i.e., the reaction order corresponds to the sum of the stoichiometric coefficients of the species appearing in the reactant of the `i` reaction. 
For example, in the reaction 20, using the matrix A defined in Eq 27, 

$$
A=\begin{bmatrix}
 1& 1 & 0 & 0  \\  \tag{51}
 0 & 0 &1  &2  
\end{bmatrix} ; 
$$

The forward reaction `i=1` has an order or 2, the backward reaction `i=2` has an order of 3. The units of the kinetic constant `k_i` of each reaction `i` can be written in terms of the order of the reaction:

$$[k_i]=\frac{M^{1-\sum_{j=1}^{s} A_{ij}}}{s} \tag{52}$$


```julia
A=([1,1,0,0],[2,1,0,0])
println("The order of the 1st reaction is $(sum(A[1]))")
println("The order of the 2nd reaction is $(sum(A[2]))")

println("The Units of k1 are s^-1 . k^$(1-(sum(A[1])))")
println("The Units of k2 are s^-1 . k^$(1-(sum(A[2])))")
```

    The order of the 1st reaction is 2
    The order of the 2nd reaction is 3
    The Units of k1 are s^-1 . k^-1
    The Units of k2 are s^-1 . k^-2


### Exercise: 
White a computer program that consists on a function that takes as input the matrices A and B for a given interaction scheme and returns: the value $(B-A)'$, the order of each reaction and the units of the rate constants. 

## 5. Differential Form of Mass Action

Now that we now how to obtain the stoichiometric matrices and the order of any given reaction, the next step is to write the differential equations.  We will use now the same formation for general interaction of eq.5 but in vector notation:

$$\begin{align*}
X_1 + X_2 \overset{k}{\longrightarrow} X_3  \tag{53}\\ 
\end{align*}$$
The stoichometrix matrix is: 


```julia
A=([1,1,0])
println("The order of the 1st reaction is $(sum(A))")
println("The Units of k are s^-1 . k^$(1-(sum(A)))")
```

    The order of the 1st reaction is 2
    The Units of k are s^-1 . k^-1


We will write down the differential equations that account for the dynamics of each species. 
$$
\begin{align}
\frac{\mathrm{d} X_1}{\mathrm{d} t}  &=&  \tag{54}\\
\frac{\mathrm{d} X_2}{\mathrm{d} t} &=& \tag{55}\\
\frac{\mathrm{d} X_3}{\mathrm{d} t} &=&  \tag{56}
\end{align}
$$

The Law of Mas Action states that the rate of change of any reactant is proportional to the concetrations of the reactants involved. The velocity depends on the affinity of the reaction and the concentrations of each of the reactants (what is before the arrow). The affinity of the reactions is defined by the rate constant `k` of the reaction. The `k` is the proportionallyty constant that correlates the speed fo the reaction with the concentrations of the species involved. In this case is $X_1$ and $X_2$:

$$
\begin{align}
\frac{\mathrm{d} X_1}{\mathrm{d} t}  &= - k X_1 X_2 \tag{57}\\
\frac{\mathrm{d} X_2}{\mathrm{d} t} &=- k X_1 X_2 \tag{58}\\
\frac{\mathrm{d} X_3}{\mathrm{d} t} &=  k X_1 X_2\tag{59} 
\end{align}
$$

Lets see if the units are correct:

$$
\begin{align}
\frac{[M]}{[s]}  &= - \frac{1}{[s][M]}[M] [M] \tag{60}\\
\frac{[M]}{[s]}  &= - \frac{1}{[s][M]}[M] [M] \tag{61}\\
\frac{[M]}{[s]}  &=  \frac{1}{[s][M]}[M] [M] \tag{62}\\
\end{align}
$$

For the case where more than one molecule of reactant is involved, we will use the reaction of formation of water. 


```julia
A=[2 1 0]
B=[0 0 2]
stoichiometric_matrix= (B-A)'
println("The order of the reaction is $(sum(A))")
println("The units of k are s^-1 . k^$(1-(sum(A)))")
println("The stoichiometric matrix is $stoichiometric_matrix ")
```

    The order of the reaction is 3
    The units of k are s^-1 . k^-2
    The stoichiometric matrix is [-2; -1; 2] 


we can decompose the equation as
 $$  X_1 +  X_1 + X_2 \overset{k}{\longrightarrow}  X_3 + X_3 \tag{63}$$
We will write down the differential equations that account for the dynamics of each species. The speed of each reaction depends on the concentration of the species involved (what is before the arrow). In this case is two units of $X_1$ and one unit of $X_2$. In addition, since the reaction consumes 2 units of $X_1$ for each unit of $X_2$, the rate of consumption of $X_1$ is twice the rate of consumption $X_2$. Another more automaticic way of doing it is to use teh matrix of stochoimetric coefficients $(B-A)'$. Therefore:
$$
\begin{align}
\frac{\mathrm{d} X_1}{\mathrm{d} t}  &= - 2k X_1^2 X_2 \tag{64}\\
\frac{\mathrm{d} X_2}{\mathrm{d} t} &=- k X_1^2 X_2 \tag{65}\\
\frac{\mathrm{d} X_3}{\mathrm{d} t} &=2  k X_1^2 X_2 \tag{66}
\end{align}
$$
Lets check if units are correct:
$$
\begin{align}
\frac{[M]}{[s]}  &=& - \frac{1}{[s][M]^2} [M]^2 [M] \tag{67}\\
\frac{[M]}{[s]}  &=& - \frac{1}{[s][M]^2} [M]^2 [M] \tag{68}\\
\frac{[M]}{[s]}  &=& \frac{1}{[s][M]^2} [M]^2 [M] \tag{69}
\end{align}
$$


## 6. General formulation of differential form of Mass action

When we have to write the differenetial equations of more complex systems of interactions, such as Eq 21, we can use a general method based omn the following equation: 
$$
\begin{align}
\frac{\mathrm{d} X}{\mathrm{d} t}&= (B-A)^T \cdot K \cdot X^A \tag{70}\\
\end{align}
$$
where  `A` and `B` are the matrices with the stoichiometric coefficients, `X` is the state vector $X=[X_1,X_2,X_3]^T$, and `K` is a matrix in the form:

$$
K=\begin{pmatrix}
 k_1 & 0   & 0  & \dots & 0 \\ 
 0 &  k_2  & 0  & \dots & 0\\ 
 0 & 0  & k_3 &\dots  &0 \\ 
 \vdots &  \vdots   &  \vdots  &  \vdots  &  \vdots  \tag{71}\\ 
 0 & 0  & 0  &  \dots & k_r 
\end{pmatrix}
$$

The calculation of $X^A$ is perfomed as:

$$X^A=\begin{pmatrix}
X_1^{A_{11}} \cdot X_2^{A_{12}} \cdot X_3^{A_{13}} \cdot  ... \cdot  X_s^{A_{1s}} \\ 
X_1^{A_{21}} \cdot X_2^{A_{22}} \cdot X_3^{A_{23}} \cdot  ... \cdot  X_s^{A_{2s}} \\ 
 \vdots  \,\,\,\,\,\,\,\,\,\,\,  \vdots  \,\,\,\,\,\,\,\,\,\,\,\,    \vdots \,\,\,\,\,\,\,\,\,   \vdots  \,\,\,\,\,\,\,\,   \vdots  \\ 
X_1^{A_{r1}} \cdot X_2^{A_{r2}} \cdot X_3^{A_{r3}} \cdot  ... \cdot  X_s^{A_{rs}} \\ 
\end{pmatrix} =\begin{pmatrix}
 \prod_{i=1}^{s} X_{i}^{A_{1i}}\\ 
\prod_{i=1}^{s} X_{i}^{A_{2i}}\\ 
   \vdots  \tag{72} \\ 
\prod_{i=1}^{s} X_{i}^{A_{ri}}\\ 
\end{pmatrix}
$$

We will use as example the reaction in eq 21. This example has two reactions $r=2$ and four species $s=4$, with the stoichiometric matrices `A` and `B` listed in Eq. 30



```julia
A=[1 1 0 0; 0 0 1 2]
B=[0 0 1 2; 1 1 0 0]
stoichiometric_matrix= (B-A)'
println("The stoichiometric matrix is $stoichiometric_matrix ")
```

    The stoichiometric matrix is [-1 1; -1 1; 1 -1; 2 -2] 


in our case

$$
K=\begin{pmatrix}
 k_1 & 0   \tag{73}\\ 
 0 &  k_2  \\ 
\end{pmatrix}
$$

and 

$$X^A=\begin{pmatrix}
X_1\cdot X_2  \tag{74}\\ 
X_3\cdot X_4^2  
\end{pmatrix}
$$

therefore

$$
\begin{align}
 \begin{bmatrix}
\frac{\mathrm{d} X_1}{\mathrm{d} t}\\ \frac{\mathrm{d} X_2}{\mathrm{d} t} \\ \frac{\mathrm{d} X_3}{\mathrm{d} t} \\ \frac{\mathrm{d} X_4}{\mathrm{d} t}
\end{bmatrix}&= \begin{bmatrix}
 -1   &  1 \\
    -1  &   1\\
     1  &  -1\\
     2  &  -2
\end{bmatrix}\begin{bmatrix}
k_1 & 0\\ 
 0& k_2
\end{bmatrix}\begin{bmatrix}
X_1 \cdot X_2\\ 
X_3 \cdot X_4^2
\end{bmatrix} \tag{75}
\end{align}
$$

so calculating

$$
\begin{align}
\begin{bmatrix}
\frac{\mathrm{d} X_1}{\mathrm{d} t}\\ \frac{\mathrm{d} X_2}{\mathrm{d} t} \\ \frac{\mathrm{d} X_3}{\mathrm{d} t} \\ \frac{\mathrm{d} X_4}{\mathrm{d} t}
\end{bmatrix}&= \begin{bmatrix}
 -1   &  1 \\
    -1  &   1\\
     1  &  -1\\
     2  &  -2
\end{bmatrix}\begin{bmatrix}
k_1 \cdot  X_1 \cdot X_2\\ 
k_2 \cdot X_3 \cdot X_4^2
\end{bmatrix} \tag{76}
\end{align}
$$

for the first equation


$$
\begin{align}
\begin{bmatrix}
\frac{\mathrm{d} X_1}{\mathrm{d} t}\\ \frac{\mathrm{d} X_2}{\mathrm{d} t} \\ \frac{\mathrm{d} X_3}{\mathrm{d} t} \\ \frac{\mathrm{d} X_4}{\mathrm{d} t}
\end{bmatrix}&= \begin{bmatrix}
 \mathbf{-1}   &  1 \\
    -1  &   1\\
     1  &  -1\\
     2  &  -2
\end{bmatrix}\begin{bmatrix}
\mathbf{k_1 \cdot  X_1 \cdot X_2}\\ 
k_2 \cdot X_3 \cdot X_4^2
\end{bmatrix} \tag{77}
\end{align}
$$

so, it is basically the Stoichiometric matrices multiplied by the vector of speeds of the reactions involved. This way, the final set of equations is:

$$\begin{align}        
            \frac{ \mathrm{d}X_1 }{\mathrm{d}t} &= - k_1 \cdot X_1 \cdot X_2 + k_2 \cdot X_3 \cdot X_4^2  \tag{78}\\ 
            \frac{ \mathrm{d}X_2 }{\mathrm{d}t} &= - k_1 \cdot X_1 \cdot X_2 + k_2 \cdot X_3 \cdot X_4^2   \tag{79}  \\
             \frac{ \mathrm{d} X_3 }{\mathrm{d}t} &=  k_1 \cdot X_1 \cdot X_2 - k_2 \cdot X_3 \cdot X_4^2  \tag{80}  \\
              \frac{ \mathrm{d} X_4 }{\mathrm{d}t} &= 2 k_1 \cdot X_1  \cdot X_2 - 2 k_2 \cdot X_3 \cdot X_4^2   \tag{81}  \\
\end{align}
$$

## 7. Mass Conservation Law

The law of conservation of mass or principle of mass conservation states that for any system closed to all transfers of matter and energy, the mass of the system must remain constant over time, as the system's mass cannot change, so quantity can neither be added nor be removed. Therefore, the quantity of mass is conserved over time. For example, in chemical reactions, the mass of the chemical components before the reaction is equal to the mass of the components after the reaction. Thus, during any chemical reaction and low-energy thermodynamic processes in an isolated system, the total mass of the reactants, or starting materials, must be equal to the mass of the products. 

We can use this law to find relationships between species of a given reaction, in such a way that we can reduce the number of variables of the systems. The mathematical formulation of thsi law is the following:

Suppose we have a solution of the system of equations defined by the stoichiometric matrix in the form $C \cdot (B-A)^T =0$ such that $C$ is a vector with number of rows $j= 1...s$ equal to the number of species involved `s`, then the following is true:

$$
\begin{align}
 C_1 \cdot X_1 + C_2 \cdot X_2 + ... + C_s \cdot X_s = cte \tag{82}
\end{align}
$$

this equation is a mass conservation law of the reaction. 

$$
\begin{align}
 C_1 \cdot X_1 + C_2 \cdot X_2 + ... + C_s \cdot X_s = cte \tag{83}
\end{align}
$$

this is true because, if we take the derivative of the last equation:
$$
\begin{align}
\frac{\mathrm{d} [C_1 \cdot X_1 + C_2 \cdot X_2 + ... + C_s \cdot X_s ]}{\mathrm{d} t}&= \frac{\mathrm{d} [C_1] \cdot X_1}{\mathrm{d} t} + \frac{\mathrm{d} [C_2] \cdot X_2}{\mathrm{d} t} +...+ \frac{\mathrm{d} [C_s] \cdot X_s}{\mathrm{d} t}  \tag{84}
\end{align}
$$
since $C$ is a vector of constant values we can write:
$$
\begin{align}
[C_1] \frac{\mathrm{d}   X_1}{\mathrm{d} t} + [C_2] \frac{\mathrm{d} X_2}{\mathrm{d} t} +...+ [C_s] \frac{\mathrm{d}  X_s}{\mathrm{d} t} = \begin{bmatrix}
C_1 & C_2 & ... & C_s \end{bmatrix} \begin{bmatrix}
\frac{\mathrm{d} X_1}{\mathrm{d} t}\\ \frac{\mathrm{d} X_2}{\mathrm{d} t} \\ \frac{\mathrm{d} X_3}{\mathrm{d} t} \\ \vdots  \\ \frac{\mathrm{d} X_s}{\mathrm{d} t} \tag{85}
\end{bmatrix}
\end{align}
$$
Now we know that, based on the general equation of the differentiatl form of mass action 

$$
\begin{align}
 \begin{bmatrix}
C_1 & C_2 & ... & C_s \end{bmatrix} \begin{bmatrix}
\frac{\mathrm{d} X_1}{\mathrm{d} t}\\ \frac{\mathrm{d} X_2}{\mathrm{d} t} \\ \frac{\mathrm{d} X_3}{\mathrm{d} t} \\ \vdots  \\ \frac{\mathrm{d} X_s}{\mathrm{d} t} \end{bmatrix} =  \begin{bmatrix}
C_1 & C_2 & ... & C_s \end{bmatrix}  (B-A)^T \cdot K \cdot X^A \tag{86}  
\end{align}
$$ 

we also know that  $C \cdot (B-A)^T =0$ by definition (vector C was a solution of the system of equations), so

$$
\begin{align}
 \begin{bmatrix}
C_1 & C_2 & ... & C_s \end{bmatrix} \begin{bmatrix}
\frac{\mathrm{d} X_1}{\mathrm{d} t}\\ \frac{\mathrm{d} X_2}{\mathrm{d} t} \\ \frac{\mathrm{d} X_3}{\mathrm{d} t} \\ \vdots  \\ \frac{\mathrm{d} X_s}{\mathrm{d} t} \end{bmatrix} =  0 \cdot K \cdot X^A \tag{87}  = 0 
\end{align}
$$ 

therefore 
$$
\begin{align}
 C_1 \cdot X_1 + C_2 \cdot X_2 + ... + C_s \cdot X_s = cte \tag{88}
\end{align}
$$


## 8. Example of Mass Conservation Law

Lets now apply the theory to find conservation laws in the reaction of the example.

$$NaCO_3 + CaCl_2 \overset{k_1}{\underset{k_2}{\longleftrightarrow}} CaCO_3 + 2 \cdot NaCl \tag{89}\\  $$

we have to find  $C \cdot (B-A)^T =0$ with $C$ being a vector with $j= 1...s$, being $s=4$. 

$$
\begin{bmatrix}
 C_1   &  C_2  &  C_3  &  C_4 \end{bmatrix}\begin{bmatrix}
 -1   &  1 \\
    -1  &   1\\
     1  &  -1\\
     2  &  -2
\end{bmatrix}=0 \tag{90}
$$

So we have to solve this system of equations:

$$\begin{cases} -C_1 - C_2 + C_3 + 2 \cdot C_4 = 0 \\ C_1 + C_2 - C_3 - 2 \cdot C_4 = 0 \end{cases} \tag{91}$$

since the two reactions correspond to a reversible system, the second equation is simply the first one multiplied by $-1$. Therefore the solytion of the system is simply
$$
C_3 + 2 \cdot C_4 - C_2 = C_1 \tag{92}
$$

so, we obtain one equation and three unknown variables. If we give values for three of the constants ($C_{1 \dots 3}=1$) and calculate the other one, we obtain:

$$
\begin{align}
1 + 2 \cdot C_4 - 1 &= 1 \tag{93}\\
C_4&=\frac{1}{2}\tag{94}
\end{align}
$$

the following the equation for the mass conservation, we obtain:


$$
\begin{align}
 X_1 + X_2 + X_3 +\frac{X_4}{2}= cte \tag{95}
\end{align}
$$

Let's test if this is correct by solving numerically the system (we will learn later how to do this). 



```julia
k1=2.3e0;  # units 1/(Ms)
k2=2.5e0;  # units 1/(M M s)

p = (k1,k2)

a₀=0.02; # units (M)
b₀=0.01; # units (M)
c₀=0; # units (M)
d₀=0; # units (M)

u₀=[a₀,b₀,c₀,d₀];
```


```julia
tspan = (0.0,300.0)
prob = ODEProblem(simpleODE3!,u₀,tspan,p)
sol = solve(prob)
plot(sol,title= "Solution of differential equations",xlabel="Time [s]",ylabel="Concentration [M]")
plot!(sol.t,sol[1,:].+sol[2,:].+sol[3,:].+0.5.*sol[4,:],label="conservation")
```




    

    



## 9. Conclusions

The law of mass action allows us to calculate the dynamcis of a set of interacting entitites of a system. It is the most widely used tool to study biologivcal systems from a quantitative predictive perspective. When possible, the mass conservation allows us to reduce the complexity of the system by reducing the number of independent variables, by providing us with dependencies between the variables of the system. This tool allows us to reduce the number of equations to solve. 
This search for mass conservation laws becomes more useful once whe have systems of more that one reversible reactions. 



## 10. Practical Exercices

-  Find the order of the reactions, the units of the kinetic rate constants, the mass conservation, and write the differential equations that govern the dynamcis of the following system:

$$ X_1  \overset{k_1}{\underset{k_2}{\longleftrightarrow}} X_2 \tag{96} $$

-  Find the order of the reactions, the units of the kinetic rate constants, the mass conservation, and write the differential equations that govern the dynamcis of the following system:
$$X_1\overset{k_1}{\longrightarrow}X_2\overset{k_2}{\longrightarrow}X_3 \tag{97}$$

-  Find the order of the reactions, the units of the kinetic rate constants, the mass conservation, and write the differential equations that govern the dynamcis of the following system:
$$
\begin{align}
 X_1 &\overset{k_1}{\longrightarrow} 2X_1  \tag{98} \\
 X_1 + X_2 &\overset{k_2}{\longrightarrow} 2X_2 \tag{99}\\
 X_2 &\overset{k_3}{\longrightarrow} 0 \tag{100}\\
 \end{align}
 $$

18ed4641d612a8d550d1415eb778fca54d1e86c0
ipynb
Jupyter Notebook
1_MassActionLaw.ipynb
davidgmiguez/julia_notebooks
b395fac8f73bf8d9d366d6354a561c722f37ce66
BSD-3-Clause
1_MassActionLaw.ipynb
davidgmiguez/julia_notebooks
b395fac8f73bf8d9d366d6354a561c722f37ce66
BSD-3-Clause
1_MassActionLaw.ipynb
davidgmiguez/julia_notebooks
b395fac8f73bf8d9d366d6354a561c722f37ce66
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES