Solving a Petersen Matrix with Julia

New to Julia
1

Julia has been great so far for my work in environmental engineering and modeling. I have only used it for simple systems thus far, but it is very performant.

Maybe I am missing something simple, but I have not found a good succinct way to input system rate equations in the form of a Petersen Matrix:
https://en.wikipedia.org/wiki/Petersen_matrix

As in the example from Wikipedia, I would normally write out the system of equations and solve

As the system gets more complex I would like to solve in matrix format:
image

But I am not sure how to do the syntax. Thanks!

2

Welcome Kester, I hope you have a great time in discourse.

Well to solve the equation you would need to use some package like DifferentialEquations.jl which I haven’t use myself.

But, lets say that you a Julia function named solve which solves the equation

\frac{\text{d}}{\text{d}t}\vec{x} = f(\vec{x}),

then in your case you have that

\vec{x} = ([A], [B], [S], [E], [ES], [P])^T.

Enough maths, lets get into Julia. The function in the right hand of the ODE, which I call “f”, has a matrix form in your case. You can write a Julia function f which takes an array x as above but also an array of parameters kgiven by

\vec{k} = (k_1, k_r, k_f, k_\mathrm{cat})^T.

Then your Julia code would look something like

f(x, k) = [-1 0 0 0;
           -2 0 0 0;
           1 -1 1 0;
           0 -1 1 1;
           0 1 -1 -1;
           0 0 0 1] * [k[1]*x[1]*x[2]^2,
                         k[2]*x[4]*x[3],
                         k[3]*x[5],
                         k[4]*x[6]]

I’m explaining myself? I hope that this is useful…

Edit: syntax correction

3

Yes that is very helpful, thank you. I have experience with simpler equations in DifferentialEquations.jl, I don’t know why I was not understanding how to do this… your example makes perfect sense though.

Part of my error was with syntax. I thought matrix format as like:

A  = [1. 0  0 -5
      4 -2  4 -3
     -4  0  0  1
      5 -2  2  3]

But I see you are using an array? Or an array of vectors? (Sorry I am still very new with data types in Julia).

EDIT: I tried the format in your answer and it gave me a dimension mismatch when attempting to use DiffEq.jl. However this format worked:

f(x, k) = [-1 0 0 0;
              -2 0 0 0; 
              1 -1 1 0; 
              0 -1 1 1; 
              0 1 -1 -1; 
              0 0 0 1] * [k[1]*x[1]*x[2]^2,
                            k[2]*x[4]*x[3],
                            k[3]*x[5],
                            k[4]*x[6]]

Still think there is more I need to learn about data types and structures in Julia. Thanks again for the help. If anyone has any more pointers it would be great!

4

Sorry, for the wrong syntax; I didn’t run the example. :sweat_smile:

5

@Kester
You can check out how to directly write reactions as shown in Wikipedia page you have cited, using Catalyst.jl package Home · Catalyst.jl and, convert it into ODESystem automatically.

6

Youre example was very helpful though, thank you!

Do you know what the difference is between the formats? I think one is a matrix, while one is a vector of vectors?

7

Thanks for the suggestion! I looked at this package, as well as ModelingToolkit.jl. I see how these would function, but I don’t quite see the benefit compared to just explicitly writing the ODEs myself and solving them. Is the goal with these packages simplicity when you have very complex systems?

With Catalyst the benefit may be based on how the problems are formulated? Since my field of research I am normally constructing a Petersen Matrix, I only need to transpose to get the system rate equation.

8

Yes you are right. You wrote a 6x4 matrix while I wrote a 6 element array of 1x4 matrices.

9

Catalyst.jl benefits for automatic conversion of reaction system to ODE, SDEs etc. If you know explicitly the matrix form it’s nice(as in your case)

But if you do not know explicitly ODEs from chemical reaction (say for large networks) , you might find Catalyst.jl beneficial as you do not have to do tedious task of checking whether the entries in matrix are correct or not. It gives you a way to explicitly find this matrix and rate expressions. Do take a look at its documentation.

1 Like
10

Simple working example

Load packages:

using DifferentialEquations
using Parameters
using StaticArrays
using Plots

Describe the diff. eq. with StaticArrays (Assume that u is a SVectorof size (6,), M a SMatrix of size (6, 4), and p a named tuple (or struct) of the all model parameters):

function f(u, p, t)
    @unpack M, k1, kf, kr, kcat = p
    A, B, S, E, ES, P = u
    v = SVector(k1*A*B^2, kf*E*S, kr*ES, kcat*ES)
    M*v
end

The reason to use StaticArrays (SVector, SMatrix, etc.) is to improve performance.

cf. https://diffeq.sciml.ai/stable/tutorials/ode_example/#ode_other_types

Parameters:

M = [
    -1  0  0  0
    -2  0  0  0
    +1 -1 +1  0
     0 -1 +1 +1
     0 +1 -1 -1
     0  0  0 +1
]

p = (
    M = SMatrix{6, 4, Float64}(M),
    k1 = 1.0, 
    kf = 1.0, 
    kr = 1.0, 
    kcat = 1.0, 
)

u0 = SVector(
    #= A  =# 1.5, 
    #= B  =# 1.2, 
    #= S  =# 0.0, 
    #= E  =# 1.0, 
    #= ES =# 0.0, 
    #= P  =# 0.0,
)

tspan = (0.0, 10.0)

Solve the problem and plot the result:

prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob)
label = ["A" "B" "S" "E" "ES" "P"]
linestyle = [:dash :dash :dashdot :dashdot :solid :solid]
plot(sol; label, linestyle, lw=1.5)

2021-09-08
2021-09-08759Ă—491 20.2 KB

You can obtain the same result with

function g(u, p, t)
    @unpack M, k1, kf, kr, kcat = p
    A, B, S, E, ES, P = u
    dA = -k1*A*B^2
    dB = -2k1*A*B^2
    dS = k1*A*B^2 - kf*E*S + kr*ES
    dE  = -kf*E*S + kr*ES + kcat*ES
    dES =  kf*E*S - kr*ES - kcat*ES
    dP = kcat*ES
    SVector(dA, dB, dS, dE, dES, dP)
end

prob = ODEProblem(g, u0, tspan, p)
sol = solve(prob)
label = ["A" "B" "S" "E" "ES" "P"]
linestyle = [:dash :dash :dashdot :dashdot :solid :solid]
plot(sol; label, linestyle, lw=1.5)

Jupyter notebook: https://github.com/genkuroki/public/blob/main/0019/Petersen%20matrix.ipynb

3 Likes
11

Thank you very much for the detailed examples, and for the info about StaticArrays. I learned a lot looking through your code here!