-
-
Notifications
You must be signed in to change notification settings - Fork 9
Description
Is your feature request related to a problem? Please describe.
I've been frustrated trying to import reaction networks in a matrix form. The documentation to load a reaction network from a matrix includes an example of the matrix hard-coded into the file, rather than a more general form that allows users to supply a .csv defining the matrix.
Describe the solution you’d like
It would be ideal to be able to supply a .csv file directly that can be converted to an appropriate matrix. I've written a minimal function that I think acomplishes this:
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{Int64}`: A transposed matrix with the stoichiometric coefficients of the substrates in each reaction.
- `rate_constants::Vector{Float64}`: A numeric vector containing the mass action 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
`julia
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])
netstoich_matrix = Matrix{Int64}(data[:, 2:end])
# Get the number of reactions (rows) and species (columns)
m_reactions, n_species = size(netstoich_matrix)
product_matrix = max.(netstoich_matrix, 0)
substrate_matrix = -min.(netstoich_matrix, 0)
# Transpose the product and substrate matrices for compatibility with catalyst.jl
product_matrix = Matrix{Int64}(transpose(product_matrix))
substrate_matrix = Matrix{Int64}(transpose(substrate_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
endThe function expects a .csv that is formatted like so:
| k | S1 | S2 | S3 |
|---|---|---|---|
| 0.1 | 1 | 0 | -1 |
| 0.25 | 0 | 1 | 1 |
| 0.12 | -1 | -1 | 0 |
With the first column k as the mass action rate constants, the column headers S1, S2.. as the species names and the matrix being the net stoichiometry matrix (only integer values allowed). This makes the structure of this .csv file very similar to writing the list of chemical reactions themselves, making this reasonably interpretable.
A possible extension could allow users to define strings in the 'k' column that would be converted to symbolic variables that could be defined later, after the reaction network is imported.
The usage would look like:
(prod, sub, k, species, m_reactions, n_species) = read_CRN("myCRN.csv")
@variables t
@species S(t)[1:n_species]
@parameters k[1:m_reactions]
# 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.rnDescribe alternatives you’ve considered
Reading the matrix from a .csv file allows users to make arbitrarily large and portable networks fairly easily. These could also be generated procedurally through some sort of loop, or by hard-coding the matrices as in the documented example.
Additional context
I had a pretty difficult time ensuring that all of the internal vectors and matrices were of the correct type to construct a ReactionSystem, so I'd like for an example like this (or even an integrated function) to be supported in the documentation for others. I am using Catalyst 15.0.8 and ReactionNetworkImporters 0.16.1
Add any other context or screenshots about the feature request here.