# Solution of a first- and second-order consecutive reaction scheme

I have a compartmental model system X, a simplified version of which, looks like this:

$$\text{System }x: \enspace \text{STO} \xrightleftharpoons[k_2]{k_1} \text{blood} \xrightarrow[]{k_3} \text{urine} \equiv x_1 \xrightleftharpoons[k_2]{k_1} x_2 \xrightarrow[]{k_3} x_3$$

The system of differential equations is then: $$\frac{dx_1}{dt} = -k_1 x_1 + k_2 x_2 \\ \frac{dx_2}{dt} = k_1 x_1 - k_2 x_2 - k_3 x_2\\ \frac{dx_3}{dt} = k_3 x_2 \\$$

I can solve this system of equations using the following matrix approach.

First, I write the rate matrix [R]. From [R] one can obtain a new matrix [A] by first replacing each diagonal element of [R] by the negative of the sum of each of row elements, and then transposing it: $$[R] = \begin{bmatrix} 0 & k_2 & 0 \\ k_1 & 0 & k_3 \\ 0 & 0 & 0 \end{bmatrix} \\ [A'] = \begin{bmatrix} -k_2 & k_2 & 0 \\ k_1 & -(k_1 + k_3) & k_3 \\ 0 & 0 & 0 \end{bmatrix} \\ [A] = \begin{bmatrix} -k_2 & k_1 & 0 \\ k_1 & -(k_1 + k_3) & 0 \\ 0 & k_3 & 0 \end{bmatrix} \\$$

I can calculate the amount in each compartment by doing the following:

$$x(t) = e^{[A]t}x(0) \enspace \text{where} \enspace x(0) = \begin{bmatrix} 0 \\ 1 \\0 \end{bmatrix}$$

In python:

RMatrix = model_matrix.as_matrix()
row, col = np.diag_indices_from(RMatrix)
RMatrix[row, col] = -(RMatrix.sum(axis=1)-RMatrix[row,col])
AMatrix = RMatrix.T

def content(t):
cont = np.dot(linalg.expm(t*AMatrix), x0))


This method is working well for me.

Now, I have a little more complicated model where reactants in compartments 1 and 2 of Systems X and Y combine to get product in System Z.

$$X + Y \rightarrow Z$$, with a reaction constant of $$k_R$$.

\begin{align} &x_1 & \xrightleftharpoons[k_2]{k_1} \enspace & x_2 &\xrightarrow[]{k_3} & x_3 \\ &+ & \ & + \\ &y_1 & \xrightleftharpoons[k_5]{k_4} \enspace & y_2 &\xrightarrow[]{k_6} & y_3 \\ &\downarrow^{k_R} & \ & \downarrow^{k_R} \\ &z_1 & \xrightleftharpoons[k_8]{k_7} \enspace & z_2 &\xrightarrow[]{k_9} & z_3 \\ \end{align}

, and the corresponding system of differential equations would be:

\begin{align} \frac{dx_1}{dt} &= -k_1 x_1 + k_2 x_2 - k_R x_1 y_1 \\ \frac{dx_2}{dt} &= k_1 x_1 - k_2 x_2 - k_3 x_2 - k_R x_2 y_2\\ \frac{dx_3}{dt} &= k_3 x_2 \\ \frac{dy_1}{dt} &= -k_4 y_1 + k_5 y_2 - k_R x_1 y_1 \\ \vdots& \\ \frac{dz_z}{dt} &= -k_7 z_1 + k_8 z_2 + k_R x_1 y_1 \\ \end{align}

I am struggling with a method to solve this system of differential equations (1st and 2nd order) to calculate the amount in each compartment at a certain time t, given the initial conditions, $$k_R$$, and the transfer rates $$k_1$$, $$k_2$$, $$k_3$$, etc...

Can I solve it using the matrix method like the one above for a system of first order differential equations? What other options in Python do I have?