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?
Thanks in advance!
