Hey Julia crew!
I am trying to import a reaction network from a .csv file to use with Catalyst.jl via Reaction NetworkImporters.jl.
I have been able to successfully construct a ReactionSystem object which can print itself using latexify, but additional commands fail with MethodErrors.
Here is a minimal example:
using Pkg
Pkg.activate(".")
using Catalyst, OrdinaryDiffEq, Plots, Latexify, CSV, DataFrames, ReactionNetworkImporters, GraphMakie, CairoMakie, NetworkLayout
function read_CRN(csvfile::String)
"""
read_CRN(csvfile::String)
Reads a CSV file containing the net reaction stoichiometry (products - reactants) and returns
two DataFrames containing substrate and product stoichiometry matrices, as well as a
numeric vector of rate constants from the first column.
### Arguments
- `csvfile::String`: The path to the input CSV file. The file should have a non-numeric
first column (rate constants) and numeric data in subsequent columns, where each row represents a reaction.
### Returns
- `product_matrix::Matrix{Int64}`: A transposed matrix with the stoichiometric coefficients of the products in each reaction.
- `substrate_matrix::Matrix{Float64}`: A transposed matrix with the stoichiometric coefficients of the substrates in each reaction.
- `rate_constants::Vector{Float64}`: A numeric vector containing the rate constants from the first column of the CSV.
- `species::Vector{String}`: A vector of species names corresponding to the columns in the stoichiometry matrix.
- `m_reactions::Int`: The number of reactions.
- `n_species::Int`: The number of species.
### Example
product_matrix, substrate_matrix, rate_constants, species, m_reactions, n_species = read_CRN("path/to/file.csv")
"""
# Read the CSV file into a DataFrame
data = CSV.read(csvfile, DataFrame)
species = names(data)[2:end]
rate_constants = Vector{Float64}(data[!, 1])
data_matrix = Matrix{Float64}(data[:, 2:end])
product_matrix = max.(data_matrix, 0)
substrate_matrix = -min.(data_matrix, 0)
# Transpose the product and substrate matrices for compatibility with catalyst.jl
product_matrix = Matrix(transpose(product_matrix))
substrate_matrix = Matrix(transpose(substrate_matrix))
# Get the number of reactions (rows) and species (columns)
m_reactions = size(data_matrix, 1) # Number of reactions (rows)
n_species = size(data_matrix, 2) # Number of species (columns)
# Convert the product matrix to integers
product_matrix = convert(Matrix{Int}, product_matrix)
# Return the transposed matrices, rate constants, species names, number of reactions, and number of species
return product_matrix, substrate_matrix, rate_constants, species, m_reactions, n_species
end
(prod, sub, k, species, m_reactions, n_species) = read_CRN("test_CRN.csv")
@species S(t)[1:n_species]
@parameters k[1:m_reactions]
@variables t
# Collect to a Vector{Num} type
pars = collect(k)
species = collect(S)
# Create the MatrixNetwork with correct types
mn = MatrixNetwork(pars, sub, prod; species=species, params=pars)
prn = loadrxnetwork(mn)
rn = prn.rn
latexify(rn)
When run on a .csv file with this contents:
| rate | S1 | S2 | S3 |
|---|---|---|---|
| 1 | -2 | 1 | 0 |
| 2 | 2 | -1 | 0 |
| 3 | -1 | -1 | 1 |
| 4 | 1 | 1 | -1 |
| 5 | 3 | 0 | -3 |
The expected latex is printed (edited here to render properly):
However, the network has not been properly constructed. Using
plot_network(rn)
returns this system:
Further, any additional calls to other functions appear to fail. Here is an example:
deficiency(rn)
MethodError: no method matching exactdiv(::Float64, ::Float64)
The function `exactdiv` exists, but no method is defined for this combination of argument types.
Closest candidates are:
exactdiv(!Matched::Integer, ::Any)
@ ModelingToolkit C:\.julia\packages\ModelingToolkit\Z9mEq\src\systems\alias_elimination.jl:376
All similar functions throw the same error (i.e. <method> exists but no method for this combination of argument types), which suggests that the ReactionSystem object was not actually constructed properly.
Does anyone have any advice for this issue?
I am using Catalyst 15.0.8. and ReactionNetworkImporters 0.16.1