Equivalent network dynamics of reaction networks using complexes (reducing the network)

[yewalenikhil65]

Based on SciML#318 and explanation from the paper
"A network dynamics approach to chemical reaction networks(CRN)" https://arxiv.org/pdf/1502.02247.pdf

Define a chemical reaction network

For demonstration purpose, lets assume a neural -stem cell regulation reaction system [1].

using Catalyst, ModelingToolkit
@parameters k[1:20] t
@variables X[1:21](t)
rxs = [Reaction(k[1]/(X[1] + (1+X[4])*(1+X[5])), [X[1]], [X[4], X[5]],[1],[1,1]),
      Reaction(k[2]/(X[1] + (1+X[4])*(1+X[5])), [X[4], X[5]], [X[1]], [1,1], [1]),
      Reaction(k[3], [X[4], X[6]], [X[7]], [1, 1], [1]),
      Reaction(k[4], [X[7]], [X[4], X[6]], [1], [1, 1]),
      Reaction(k[5], [X[7]], [X[8]], [1], [1]),
      Reaction(k[6], [X[8]], [X[7]], [1], [1]),
      Reaction(k[7], [X[9] , X[10]], [X[11]], [1,1],[1]),
      Reaction(k[8], [X[11]], [X[9] , X[10]], [1],[1,1]),
      Reaction(k[9], [X[2], X[12]],[X[13]], [1,1], [1]),
      Reaction(k[10], [X[13]], [X[2], X[12]], [1], [1,1]),
      Reaction(k[11], [X[13]], [X[8]],[1], [1]),
      Reaction(k[12], [X[8]], [X[13]], [1], [1]),
      Reaction(k[13], [X[14], X[15]],[X[16]], [1,1],[1]),
      Reaction(k[14], [X[16]], [X[14], X[15]], [1],[1,1]),
      Reaction(k[15]*X[16]/(1 + X[17] + X[18] + X[19] + X[18]*X[19]), [X[17], X[18]],[X[19]], [1,1],[1]),
      Reaction(k[16]*X[16]/(1 + X[17] + X[18] + X[19] + X[18]*X[19]), [X[19]],[X[17],X[18]], [1],[1,1]),
      Reaction(k[17]/(1 + X[19]), [X[20], X[21]],[X[3]], [1,1], [1]),
      Reaction(k[18]/(1 + X[19]), [X[3]],[X[20], X[21]], [1], [1,1]),
      Reaction(k[19], [X[21]], [X[8]], [1], [1]),
      Reaction(k[20], [X[8]], [X[21]], [1], [1])
     ]

rs = ReactionSystem(rxs, t, collect(X), collect(k))

# Rate constants 
ps = [√2*5/101; √2*5/101; 1; 1; 1; 1; 1; 1;
    1/101; 1/101; 1; 1; 1; 1; 1; 1;
    1/101; 1/101; 1;1];

# and initial condition
x0 = [5.0; 5; 0; 0; 0; 5; 0; 0; 5; 5; 0;
        5; 0; 5; 5; 0; 5; 5; 0; 5; 5]

tspan = (0.0, 550.0)
ratelaw = oderatelaw.(equations(rs))

Some useful stoichiomatrices based on "complex" representation of the CRN based on [2],

complexes_in_reactions = reaction_complexes(rs); # unique complexes in the network
kd = reaction_rates(rs);   # reaction rate
KD = diagm((map(Num, kd)));  # reaction rates diagonal matrix

Z = complex_stoich_matrix(rs); # k'th column of Z indicates composition of k'th complex.

B = complex_incidence_matrix(rs); 
#=
Any CRN is
defined by a c × r incidence matrix B,
B[i,j] =   -1 ,if the ith complex is the substrate of the jth reaction,
            1 ,if the ith complex is the product of the jth reaction,
            0 ,otherwise.
=#

Δ = complex_outgoing_matrix(rs);
#=
Any CRN is defined by a c × r Outgoing matrix B,
Δ[i,j] =   0 ,if B[i,j] = 1
         B[i,j] , otherwise
=#

S = Z*B; # net stoichiometric matrix

ns = numspecies(rs);   # number of species in reactionSystem     
nr = numreactions(rs);  # number of reactions in reactionSystem
nc = size(complexes_in_reactions,1); # number of complexes in the network

Visualise original complex CRN as follows

# original network as complexes
gr()
plot_network(complexes_in_reactions, B)

As stated in [2], CRN can be re-wrttin in ODE form as . Thus, we can express CRN Network as ODESystem as follows.

expr = build_function(KD, [species(rs); parameters(rs)])
out = zeros(Float64, nr, nr);
f = eval(expr[2]);

# original model or oprob = convert(ODEProblem, rs)
function dudt!(du,u,p,t)
    f(out, [u;ps])
    du .= -Z*B*out*Δ'*exp.(Z'*log.(u))
end

prob = ODEProblem(dudt!, x0 .+ eps(), tspan, ps)  # eps() for avoiding log.(0)
sol = solve(prob)
Sol = Array(sol)

Following procedure outputs the index of species to be eliminated from the network.

• Extract index set corresponding to complexes in the network from complex composition matrix Z.
• For given pair of complexes and that share atleast one species ,define = () \ {α}. If they don't share any species set = 𝛷
• Define = as indices to be eliminated.

## steps to reduce the reactions
function SpeciesToEliminate(Z::Array{Int64, 2})
    Index_set = []
    for i ∈ 1:nc
        push!(Index_set, convert(Vector{Int64},
                    indexin(complexes_in_reactions[:,1][i][:,1], X)))
    end
    Eliminate_species = Int64[]
    no_species = Int64[]
    for i ∈ 1:nc-1
        for j ∈ i+1:nc
            Complex_pair = [i, j];
            Shared_species = intersect(Index_set[i], Index_set[j]);
            if  length(Shared_species) != 0
                El₁ = setdiff(Index_set[i], Shared_species);
                El₂ = setdiff(Index_set[j], Shared_species);
                El = union(El₁, El₂);
                append!(Eliminate_species, union(Eliminate_species, El));
            end
            Complex_rest = setdiff(1:nc, Complex_pair);
            for n ∈ 1:length(Complex_rest)
                m = Complex_rest[n];
                if length(setdiff(Index_set[m], Eliminate_species)) == 0
                   append!(no_species ,union(no_species, Index_set[m]));
                end
            end
        end
    end
    Eliminate_species = setdiff(Eliminate_species, no_species);
end
I = SpeciesToEliminate(Z)

Following procedure eliminates species identified from previous step

• Replace rows of net-stoichiomatrix corresponding to the elements \ β with zeros and call this Sᵦ, for every β ∈
• Find ker(Sᵦᵀ), and use appropriate columns of these matrix to join certain reactions into a single linkage class according to the Preposition 1. in [1]

function Elimination_of_species(x::Vector{Float64})
    I = SpeciesToEliminate(Z)
    Elimination = []
    for j ∈ 1:length(I)
        p = I[j];
        I₁ = setdiff(I, p);
        n₁ = length(I₁);
        S₁ = copy(S);
        S₁[I₁,:] = zeros(n₁, nr);
        N = rational_basis_nullspace(S₁');
        cl = size(N,2);
        M = Matrix(undef, size(N,1), 0)
        for i ∈ 1:cl
            if N[p,i] != 0
                M = hcat(M, N[:,i])
                break
            end
        end
        C = dot(M, x0);
        a = M[p];
        M[p] = 0;
        push!(Elimination, (-dot(M,x) + C)/a);
    end
    return Elimination
end

Certain complexes then become identical. So we delete all the identical columns from Zₑ and keep only one of them.

Deleted_identical_complexes = [];
Remaining_identical_complexes = [];
for i = 1:nc-1
    for j = i+1:nc
        if (Zₑ[:,i] == Zₑ[:,j] && maximum(Zₑ[:,j]) != 0)
             push!(Deleted_identical_complexes, j);
             push!(Remaining_identical_complexes, i);
        end
    end
end
Zₑ = Zₑ[:, setdiff(1:end, Deleted_identical_complexes)]

We determine the ncₑ x nrₑ incidence composition matrix Bₑ of the equivalent network.

Bₑ = copy(B);
y = length(Remaining_identical_complexes);
z = length(Deleted_identical_complexes);
for i ∈ 1:z
    m = Deleted_identical_complexes[i];
    p = Remaining_identical_complexes[i];
    for j ∈ 1:nr
        if (Bₑ[p,j] == 0 && Bₑ[m,j] != 0)
           Bₑ[p,j] = Bₑ[m,j];
        end
    end
end
Bₑ = Bₑ[setdiff(1:end, Deleted_identical_complexes), :];
nsₑ, ncₑ = size(Zₑ);
nrₑ = size(Bₑ,2);
Δₑ = copy(Bₑ)
Δₑ[Δₑ .== 1] .= 0  # complex outgoing matrix as defined previously

## equivalent network complexes in reactions
complexes_in_reactionsₑ = complexes_in_reactions[setdiff((1:end), Deleted_identical_complexes),:][:,1]
comp_reacₑ = []
for i ∈ 1:ncₑ
    ind = findall(in(collect(X)[I]), complexes_in_reactionsₑ[i][:,1])
    if ind != Int64[]  # didnt found the species
         complexes_in_reactionsₑ[i] = complexes_in_reactionsₑ[i][setdiff(1:end, ind), :];
    end
    push!(comp_reacₑ, sum(complexes_in_reactionsₑ[i][:,1].*complexes_in_reactionsₑ[i][:,2]))
end
complexes_in_reactionsₑ = [complexes_in_reactionsₑ comp_reacₑ]
plot_network(complexes_in_reactionsₑ, Bₑ)

Equivalent reduced model after elimination

# Equivalent /reduced model
function dudtᵣₑ!(du,u,p,t)
    z = p[1]
    b = p[2]
    δ = p[3]
    reduced = p[4]
    deleted = p[5]
    Schur,u = SchurComplement(b,δ,u,reduced,deleted)
    du .= -z[:,reduced]*Schur*exp.(z[:,reduced]'*log.(u))
end

solᵣₑ = solve(ODEProblem(dudtᵣₑ!, x0 .+ eps(), tspan,
                                        [Z,B,Δ,setdiff(1:nc, Deleted_identical_complexes), Deleted_identical_complexes ]))

Plot for species which are retained, it matches fairly well(qualitatively) with original model. Lets check for X[2] (t) species,

plot(sol.t, sol[2,:], color =:red, line = 4, label = "Original Model")
plot!(solᵣₑ.t, solᵣₑ[2,:], color =:blue, line = 4, label = "Reduced Model")
xlabel!("Time(sec)")
ylabel!("X₂(t)")

Lets plot for species which is eliminated, it doesn't match the original model as it has been eliminated from the calculation,and thus retaining its initial concentration

plot(sol.t, sol[6,:], color =:red, line = 4, label = "Original Model")
plot!(solᵣₑ.t, solᵣₑ[6,:], color =:blue, line = 4, label = "Reduced Model")
xlabel!("Time(sec)")
ylabel!("X₆(t)")

References

1. An Automated Model Reduction Method for Biochemical Reaction Networks ,Manvel Gasparyan, Arnout Van Messem ,Shodhan Rao ,https://www.mdpi.com/2073-8994/12/8/1321
2. "A network dynamics approach to chemical reaction networks(CRN)" https://arxiv.org/pdf/1502.02247.pdf

[yewalenikhil65]

Functions

reaction_complexes()
reaction_rates()
complex_stoich_matrix()
complex_incidence_matrix()
complex_outgoing_matrix()

are as defined in PR SciML#359

[isaacsas]

@yewalenikhil65 these seem like a good set of API functions. The only name I don't like is reaction_rates, which is too vague I think. I'll add specific comments on their implementations in the PR.

[isaacsas]

I haven't read the model reduction approach you cite, so can't comment on that method yet.

[yewalenikhil65]

@yewalenikhil65 these seem like a good set of API functions. The only name I don't like is reaction_rates, which is too vague I think. I'll add specific comments on their implementations in the PR.

@isaacsas
Sure, we can use any other function name as you suggest.

I also have some function to plot the network with nodes as complexes and edges as reaction arrows. If that is also good for catalyst, I would like to add them. It's the plot_network function from above script

[isaacsas]

Yeah, that sounds great. Let's add that as a second PR though after we get the first merged. (Just to keep it from getting much longer, which will make it take longer for me to review.)

[yewalenikhil65]

Sure