Skip to main content

gapsmith_core/
model.rs

1//! In-memory metabolic model.
2//!
3//! Analogue of cobrar's `ModelOrg` S4 class. Replaces RDS serialization with
4//! serde (CBOR/JSON — see `gapsmith-io`).
5
6use crate::{Compartment, CpdId, GeneId, Metabolite, Reaction, RxnId, StoichMatrix};
7use serde::{Deserialize, Serialize};
8use std::collections::HashMap;
9
10/// Provenance and version metadata that travels with the model.
11#[derive(Clone, Debug, Default, Serialize, Deserialize)]
12pub struct ModelAnnot {
13    /// Genome / organism identifier (usually the FASTA basename).
14    pub id: String,
15    /// Human-readable name.
16    #[serde(default, skip_serializing_if = "Option::is_none")]
17    pub name: Option<String>,
18    /// gapseq version that produced the model.
19    #[serde(default, skip_serializing_if = "Option::is_none")]
20    pub gapsmith_version: Option<String>,
21    /// Reference sequence database version.
22    #[serde(default, skip_serializing_if = "Option::is_none")]
23    pub seqdb_version: Option<String>,
24    /// Taxonomic domain (`Bacteria`, `Archaea`, ...).
25    #[serde(default, skip_serializing_if = "Option::is_none")]
26    pub tax_domain: Option<String>,
27    /// Gram staining (`pos`, `neg`, `na`).
28    #[serde(default, skip_serializing_if = "Option::is_none")]
29    pub gram: Option<String>,
30    /// Free-form notes (merged into SBML `<notes>` on export).
31    #[serde(default)]
32    pub notes: Vec<String>,
33}
34
35#[derive(Debug, thiserror::Error)]
36pub enum ModelError {
37    #[error("unknown reaction id `{0}`")]
38    UnknownRxn(String),
39    #[error("unknown metabolite id `{0}`")]
40    UnknownMet(String),
41    #[error("duplicate reaction id `{0}`")]
42    DuplicateRxn(String),
43    #[error("duplicate metabolite id `{0}`")]
44    DuplicateMet(String),
45    #[error(
46        "stoichiometric matrix shape {got:?} does not match model dimensions {expected:?}"
47    )]
48    ShapeMismatch { got: (usize, usize), expected: (usize, usize) },
49}
50
51#[derive(Clone, Debug, Default, Serialize, Deserialize)]
52pub struct Model {
53    pub annot: ModelAnnot,
54    pub compartments: Vec<Compartment>,
55    pub mets: Vec<Metabolite>,
56    pub rxns: Vec<Reaction>,
57    /// Genes referenced by any reaction's GPR.
58    #[serde(default)]
59    pub genes: Vec<GeneId>,
60    /// Sparse stoichiometric matrix in CSC (one column per reaction).
61    pub s: StoichMatrix,
62}
63
64impl Model {
65    /// Empty model with the three standard compartments (`c0`, `e0`, `p0`).
66    pub fn new(id: impl Into<String>) -> Self {
67        Self {
68            annot: ModelAnnot { id: id.into(), ..Default::default() },
69            compartments: Compartment::default_three(),
70            mets: Vec::new(),
71            rxns: Vec::new(),
72            genes: Vec::new(),
73            s: StoichMatrix::zeros(0, 0),
74        }
75    }
76
77    pub fn rxn_count(&self) -> usize {
78        self.rxns.len()
79    }
80    pub fn met_count(&self) -> usize {
81        self.mets.len()
82    }
83
84    /// Lookup map from reaction id → index. Built on demand (O(n)); callers
85    /// that need many lookups should cache the result.
86    pub fn rxn_index(&self) -> HashMap<RxnId, usize> {
87        self.rxns
88            .iter()
89            .enumerate()
90            .map(|(i, r)| (r.id.clone(), i))
91            .collect()
92    }
93
94    pub fn met_index(&self) -> HashMap<CpdId, usize> {
95        self.mets
96            .iter()
97            .enumerate()
98            .map(|(i, m)| (m.id.clone(), i))
99            .collect()
100    }
101
102    /// Sanity-check that the stoichiometric matrix shape matches the model.
103    pub fn check_shape(&self) -> Result<(), ModelError> {
104        let expected = (self.mets.len(), self.rxns.len());
105        let got = (self.s.rows(), self.s.cols());
106        if got == expected {
107            Ok(())
108        } else {
109            Err(ModelError::ShapeMismatch { got, expected })
110        }
111    }
112}
113
114#[cfg(test)]
115mod tests {
116    use super::*;
117    use crate::{CompartmentId, Reversibility};
118
119    fn toy_model() -> Model {
120        // 3 mets × 2 rxns: A -> B, B -> C.
121        let mut m = Model::new("toy");
122        m.mets.push(Metabolite::new("cpdA", "A", CompartmentId::CYTOSOL));
123        m.mets.push(Metabolite::new("cpdB", "B", CompartmentId::CYTOSOL));
124        m.mets.push(Metabolite::new("cpdC", "C", CompartmentId::CYTOSOL));
125
126        let mut r1 = Reaction::new("r1", "A -> B", 0.0, 1000.0);
127        r1.obj_coef = 0.0;
128        m.rxns.push(r1);
129        m.rxns.push(Reaction::new("r2", "B -> C", 0.0, 1000.0));
130
131        m.s = StoichMatrix::from_triplets(
132            3,
133            2,
134            vec![(0, 0, -1.0), (1, 0, 1.0), (1, 1, -1.0), (2, 1, 1.0)],
135        );
136        m
137    }
138
139    #[test]
140    fn shape_check() {
141        let m = toy_model();
142        m.check_shape().unwrap();
143    }
144
145    #[test]
146    fn bad_shape_detected() {
147        let mut m = toy_model();
148        m.rxns.push(Reaction::new("r3", "C -> sink", 0.0, 1000.0));
149        assert!(m.check_shape().is_err());
150    }
151
152    #[test]
153    fn reaction_indexing() {
154        let m = toy_model();
155        let idx = m.rxn_index();
156        assert_eq!(idx[&RxnId::new("r1")], 0);
157        assert_eq!(idx[&RxnId::new("r2")], 1);
158    }
159
160    #[test]
161    fn reversibility_from_bounds() {
162        let m = toy_model();
163        assert_eq!(m.rxns[0].reversibility(), Reversibility::Forward);
164    }
165}