Connected Reversible Linear Reaction Solver

What is being modeled?

This system is taken from 4.4. Connected Reversible Linear Reactions in [1]. I have rewritten the chemical equation to simplify the notation of the fluxes over the reactions.

\[ x_1 \underset{v_2}{\stackrel{v_1}{\rightleftharpoons}} x_2 \stackrel{v_3} \rightarrow x_3 \underset{v_5}{\stackrel{v_4}{\rightleftharpoons}} x_4 \]

\[ \mathbf S = \begin{pmatrix} {-1} & {1} & {0} & {0} & {0} \\ {1} & {-1} & {-1} & {0} & {0} \\ {0} & {0} & {1} & {-1} & {1} \\ {0} & {0} & {0} & {1} & {-1} \\ \end{pmatrix} \]

\[ \mathbf v = \begin{pmatrix} v_1 \\ v_2 \\ v_3 \\ v_4 \\ v_5 \end{pmatrix} \]

\[ \mathbf{Sv} = \mathbf 0 \]

Steady state fluxes

This demo will find steady state fluxes using JuMP.jl [2] with a GLPK solver.

Create the model

# Use JuMP.jl with GLPK
model = Model(GLPK.Optimizer)

# Define flux variables and their constraints
@variable(model, 0.0 <= v1 <= 1.0)
@variable(model, 0.0 <= v2 <= 1.0)
@variable(model, 0.0 <= v3 <= 1.0)
@variable(model, 0.0 <= v4 <= 1.0)
@variable(model, 0.0 <= v5 <= 1.0)

# Steady state constraints, Sv=0
@constraint(model, -v1 + v2 == 0)
@constraint(model, v1 - v2 - v3 == 0)
@constraint(model, v3 - v4 + v5 == 0)
@constraint(model, v4 - v5 == 0)

# Set objective function to maximize
@objective(model, Max, v4)

# Solve the model
optimize!(model)

Results of optimization

Print diagnostic information.

println("Was solution feasible? ", is_solved_and_feasible(model))

Was solution feasible? true

What was the solution?

println("Optimal v1: ", value(v1))
println("Optimal v2: ", value(v2))
println("Optimal v3: ", value(v3))
println("Optimal v4: ", value(v4))
println("Optimal v5: ", value(v5))
println("Optimal objective value: ", objective_value(model))

Optimal v1: 0.0
Optimal v2: 0.0
Optimal v3: 0.0
Optimal v4: 1.0
Optimal v5: 1.0
Optimal objective value: 1.0

Dynamic model

This plot reproduces Fig. 4.3 in [1]. It is the time course of concentrations given the equations below, and starting with the given initial conditions.

I am writing the chemical equation to explcitly show the rate constants used in the vector below:

\[ x_1 \underset{k_{-1}}{\stackrel{k_1}{\rightleftharpoons}} x_2 \stackrel{k_2} \rightarrow x_3 \underset{k_{-3}}{\stackrel{k_3}{\rightleftharpoons}} x_4 \]

The stoichiometric is the same as above. I am writing the flux vector differently than above to show the concentrations and rate constants explicitly.

\[ \mathbf S = \begin{pmatrix} {-1} & {1} & {0} & {0} & {0} \\ {1} & {-1} & {-1} & {0} & {0} \\ {0} & {0} & {1} & {-1} & {1} \\ {0} & {0} & {0} & {1} & {-1} \\ \end{pmatrix}, \mathbf v(\mathbf x) = \begin{pmatrix} k_1 x_1 \\ k_{-1} x_2 \\ k_2 x_2 \\ k_3 x_3 \\ k_{-3}x_4 \end{pmatrix} \]

This gives the following differential equations for concentrations over time:

\[ {dx_1 \over dt} = -k_1x_1+k_{-1}x_2 \] \[ {dx_2 \over dt} = k_1x_1-k_{-1}x_2-k_2x_2 \] \[ {dx_3 \over dt} = k_2x_2-k_3x_3+k_{-3}x_4 \] \[ {dx_4 \over dt} = k_3x_3-k_{-3}x_4 \]

Compute concentrations over time

Set initial conditions and compute concentrations over time.

# Rate constants
k_fwd_1 = 1.0
k_rev_1 = 1.0
k_fwd_2 = 1.0
k_fwd_3 = 1.0
k_rev_3 = 1.0

# Initialize matrix that will store trajectory
x = zeros(Float64, 1000, 4)

# Set initial concentration conditions
x[1, 1] = 1.0
x[1, 2] = 0.0
x[1, 3] = 0.0
x[1, 4] = 0.0

# 10 arbitrary units of time. Set x axis for plotting.
t = 10
ts = range(start=0, stop=t, length=length(x[:, 1]))

# Step size
h = t / length(x[:, 1])

# Comptue concentration trajectories
for i ∈ 2:length(x[:, 1])
    x1 = x[i-1, 1]
    x2 = x[i-1, 2]
    x3 = x[i-1, 3]
    x4 = x[i-1, 4]
    x[i, 1] = x1 + h*(-k_fwd_1*x1 + k_rev_1*x2)
    x[i, 2] = x2 + h*(k_fwd_1*x1 - k_rev_1*x2 - k_fwd_2*x2)
    x[i, 3] = x3 + h*(k_fwd_2*x2 - k_fwd_3*x3 + k_rev_3*x4)
    x[i, 4] = x4 + h*(k_fwd_3*x3 - k_rev_3*x4)
end

Plot concentrations over time

plot(
    ts,
    x, 
    labels=["x1" "x2" "x3" "x4"],
    linewidth=3, 
    xlabel="Time",
    ylabel="Concentration"
)

Citation

Bernhard Ø. Palsson. Systems Biology: Simulation of Dynamic Network States. Cambridge University Press, 2011. doi:10.1017/CBO9780511736179.
Lubin, M. et al. JuMP 1.0: recent improvements to a modeling language for mathematical optimization. Math. Prog. Comp. 15, 581–589 (2023).