oldies_copasi/
lib.rs

1//! # COPASI-RS
2//!
3//! Revival of COPASI biochemical network simulator in Rust.
4//!
5//! ## History
6//!
7//! COPASI (Complex Pathway Simulator) was developed at the Virginia
8//! Bioinformatics Institute and EML Research. It simulates biochemical
9//! networks using ODEs, stochastic methods, and hybrid approaches.
10//!
11//! ## SBML Support
12//!
13//! This crate also provides SBML (Systems Biology Markup Language) import
14//! capabilities, the standard format for biochemical models.
15//!
16//! ## Features
17//!
18//! 1. **ODE Simulation**: Deterministic simulation with LSODA
19//! 2. **Stochastic**: Gillespie's SSA (Stochastic Simulation Algorithm)
20//! 3. **Hybrid**: Adaptive switching between deterministic/stochastic
21//! 4. **Steady State**: Newton's method for equilibrium
22//! 5. **Parameter Estimation**: Levenberg-Marquardt, genetic algorithms
23//! 6. **Sensitivity Analysis**: Local and global sensitivity
24
25use oldies_core::{OldiesError, Result, Time};
26use ndarray::{Array1, Array2};
27use serde::{Deserialize, Serialize};
28use std::collections::HashMap;
29
30// =============================================================================
31// SBML CORE TYPES
32// =============================================================================
33
34/// SBML Level/Version
35#[derive(Debug, Clone, Copy, Serialize, Deserialize)]
36pub struct SbmlVersion {
37    pub level: u8,
38    pub version: u8,
39}
40
41impl Default for SbmlVersion {
42    fn default() -> Self {
43        Self { level: 3, version: 2 }
44    }
45}
46
47/// SBML Unit definition
48#[derive(Debug, Clone, Serialize, Deserialize)]
49pub struct UnitDefinition {
50    pub id: String,
51    pub name: Option<String>,
52    pub units: Vec<Unit>,
53}
54
55/// SBML Unit
56#[derive(Debug, Clone, Serialize, Deserialize)]
57pub struct Unit {
58    pub kind: UnitKind,
59    pub exponent: f64,
60    pub scale: i32,
61    pub multiplier: f64,
62}
63
64/// Standard unit kinds
65#[derive(Debug, Clone, Copy, Serialize, Deserialize)]
66pub enum UnitKind {
67    Mole,
68    Litre,
69    Second,
70    Metre,
71    Kilogram,
72    Item,
73    Dimensionless,
74}
75
76/// Compartment (reaction container)
77#[derive(Debug, Clone, Serialize, Deserialize)]
78pub struct Compartment {
79    pub id: String,
80    pub name: Option<String>,
81    pub spatial_dimensions: u8,
82    pub size: f64,
83    pub units: Option<String>,
84    pub constant: bool,
85}
86
87impl Compartment {
88    pub fn new(id: &str, size: f64) -> Self {
89        Self {
90            id: id.to_string(),
91            name: None,
92            spatial_dimensions: 3,
93            size,
94            units: None,
95            constant: true,
96        }
97    }
98}
99
100/// Species (molecule, protein, metabolite)
101#[derive(Debug, Clone, Serialize, Deserialize)]
102pub struct Species {
103    pub id: String,
104    pub name: Option<String>,
105    pub compartment: String,
106    pub initial_concentration: Option<f64>,
107    pub initial_amount: Option<f64>,
108    pub substance_units: Option<String>,
109    pub has_only_substance_units: bool,
110    pub boundary_condition: bool,
111    pub constant: bool,
112}
113
114impl Species {
115    pub fn new(id: &str, compartment: &str, initial_concentration: f64) -> Self {
116        Self {
117            id: id.to_string(),
118            name: None,
119            compartment: compartment.to_string(),
120            initial_concentration: Some(initial_concentration),
121            initial_amount: None,
122            substance_units: None,
123            has_only_substance_units: false,
124            boundary_condition: false,
125            constant: false,
126        }
127    }
128}
129
130/// Parameter (kinetic constant)
131#[derive(Debug, Clone, Serialize, Deserialize)]
132pub struct Parameter {
133    pub id: String,
134    pub name: Option<String>,
135    pub value: f64,
136    pub units: Option<String>,
137    pub constant: bool,
138}
139
140impl Parameter {
141    pub fn new(id: &str, value: f64) -> Self {
142        Self {
143            id: id.to_string(),
144            name: None,
145            value,
146            units: None,
147            constant: true,
148        }
149    }
150}
151
152/// Species reference in a reaction
153#[derive(Debug, Clone, Serialize, Deserialize)]
154pub struct SpeciesReference {
155    pub species: String,
156    pub stoichiometry: f64,
157    pub constant: bool,
158}
159
160impl SpeciesReference {
161    pub fn new(species: &str, stoichiometry: f64) -> Self {
162        Self {
163            species: species.to_string(),
164            stoichiometry,
165            constant: true,
166        }
167    }
168}
169
170/// Kinetic law expression
171#[derive(Debug, Clone, Serialize, Deserialize)]
172pub enum KineticLaw {
173    /// Mass action: k * [A]^a * [B]^b
174    MassAction {
175        rate_constant: String,
176    },
177    /// Michaelis-Menten: Vmax * [S] / (Km + [S])
178    MichaelisMenten {
179        vmax: String,
180        km: String,
181        substrate: String,
182    },
183    /// Hill equation: Vmax * [S]^n / (K^n + [S]^n)
184    Hill {
185        vmax: String,
186        k: String,
187        substrate: String,
188        n: f64,
189    },
190    /// Reversible Michaelis-Menten
191    ReversibleMM {
192        vmax_f: String,
193        km_f: String,
194        vmax_r: String,
195        km_r: String,
196    },
197    /// Custom expression (MathML string or infix)
198    Custom(String),
199}
200
201/// Reaction
202#[derive(Debug, Clone, Serialize, Deserialize)]
203pub struct Reaction {
204    pub id: String,
205    pub name: Option<String>,
206    pub reversible: bool,
207    pub reactants: Vec<SpeciesReference>,
208    pub products: Vec<SpeciesReference>,
209    pub modifiers: Vec<String>,
210    pub kinetic_law: KineticLaw,
211    pub local_parameters: Vec<Parameter>,
212}
213
214impl Reaction {
215    /// Create a simple A -> B reaction
216    pub fn simple(id: &str, reactant: &str, product: &str, rate_constant: &str) -> Self {
217        Self {
218            id: id.to_string(),
219            name: None,
220            reversible: false,
221            reactants: vec![SpeciesReference::new(reactant, 1.0)],
222            products: vec![SpeciesReference::new(product, 1.0)],
223            modifiers: Vec::new(),
224            kinetic_law: KineticLaw::MassAction {
225                rate_constant: rate_constant.to_string(),
226            },
227            local_parameters: Vec::new(),
228        }
229    }
230
231    /// Create an enzymatic reaction with Michaelis-Menten kinetics
232    pub fn enzymatic(
233        id: &str,
234        substrate: &str,
235        product: &str,
236        enzyme: &str,
237        vmax: &str,
238        km: &str,
239    ) -> Self {
240        Self {
241            id: id.to_string(),
242            name: None,
243            reversible: false,
244            reactants: vec![SpeciesReference::new(substrate, 1.0)],
245            products: vec![SpeciesReference::new(product, 1.0)],
246            modifiers: vec![enzyme.to_string()],
247            kinetic_law: KineticLaw::MichaelisMenten {
248                vmax: vmax.to_string(),
249                km: km.to_string(),
250                substrate: substrate.to_string(),
251            },
252            local_parameters: Vec::new(),
253        }
254    }
255}
256
257/// Assignment rule (algebraic constraint)
258#[derive(Debug, Clone, Serialize, Deserialize)]
259pub struct AssignmentRule {
260    pub variable: String,
261    pub expression: String,
262}
263
264/// Rate rule (ODE definition)
265#[derive(Debug, Clone, Serialize, Deserialize)]
266pub struct RateRule {
267    pub variable: String,
268    pub expression: String,
269}
270
271/// Event (discrete state change)
272#[derive(Debug, Clone, Serialize, Deserialize)]
273pub struct Event {
274    pub id: String,
275    pub trigger: String,
276    pub delay: Option<f64>,
277    pub assignments: Vec<EventAssignment>,
278}
279
280/// Event assignment
281#[derive(Debug, Clone, Serialize, Deserialize)]
282pub struct EventAssignment {
283    pub variable: String,
284    pub expression: String,
285}
286
287// =============================================================================
288// SBML MODEL
289// =============================================================================
290
291/// Complete SBML model
292#[derive(Debug, Clone, Serialize, Deserialize)]
293pub struct SbmlModel {
294    pub id: String,
295    pub name: Option<String>,
296    pub sbml_version: SbmlVersion,
297    pub compartments: Vec<Compartment>,
298    pub species: Vec<Species>,
299    pub parameters: Vec<Parameter>,
300    pub reactions: Vec<Reaction>,
301    pub assignment_rules: Vec<AssignmentRule>,
302    pub rate_rules: Vec<RateRule>,
303    pub events: Vec<Event>,
304}
305
306impl SbmlModel {
307    pub fn new(id: &str) -> Self {
308        Self {
309            id: id.to_string(),
310            name: None,
311            sbml_version: SbmlVersion::default(),
312            compartments: Vec::new(),
313            species: Vec::new(),
314            parameters: Vec::new(),
315            reactions: Vec::new(),
316            assignment_rules: Vec::new(),
317            rate_rules: Vec::new(),
318            events: Vec::new(),
319        }
320    }
321
322    /// Add a compartment
323    pub fn add_compartment(&mut self, compartment: Compartment) {
324        self.compartments.push(compartment);
325    }
326
327    /// Add a species
328    pub fn add_species(&mut self, species: Species) {
329        self.species.push(species);
330    }
331
332    /// Add a parameter
333    pub fn add_parameter(&mut self, parameter: Parameter) {
334        self.parameters.push(parameter);
335    }
336
337    /// Add a reaction
338    pub fn add_reaction(&mut self, reaction: Reaction) {
339        self.reactions.push(reaction);
340    }
341
342    /// Get species by ID
343    pub fn get_species(&self, id: &str) -> Option<&Species> {
344        self.species.iter().find(|s| s.id == id)
345    }
346
347    /// Get parameter by ID
348    pub fn get_parameter(&self, id: &str) -> Option<&Parameter> {
349        self.parameters.iter().find(|p| p.id == id)
350    }
351
352    /// Build stoichiometry matrix
353    pub fn stoichiometry_matrix(&self) -> Array2<f64> {
354        let n_species = self.species.len();
355        let n_reactions = self.reactions.len();
356        let mut matrix = Array2::zeros((n_species, n_reactions));
357
358        let species_index: HashMap<_, _> = self.species.iter()
359            .enumerate()
360            .map(|(i, s)| (s.id.clone(), i))
361            .collect();
362
363        for (j, reaction) in self.reactions.iter().enumerate() {
364            // Reactants (negative stoichiometry)
365            for sr in &reaction.reactants {
366                if let Some(&i) = species_index.get(&sr.species) {
367                    matrix[[i, j]] -= sr.stoichiometry;
368                }
369            }
370            // Products (positive stoichiometry)
371            for sr in &reaction.products {
372                if let Some(&i) = species_index.get(&sr.species) {
373                    matrix[[i, j]] += sr.stoichiometry;
374                }
375            }
376        }
377
378        matrix
379    }
380}
381
382// =============================================================================
383// SIMULATOR
384// =============================================================================
385
386/// Simulation method
387#[derive(Debug, Clone, Copy)]
388pub enum SimulationMethod {
389    /// Deterministic ODE (LSODA)
390    Deterministic,
391    /// Stochastic (Gillespie SSA)
392    Stochastic,
393    /// Hybrid (adaptive switching)
394    Hybrid,
395    /// Tau-leaping (approximate stochastic)
396    TauLeaping,
397}
398
399/// Simulation result
400#[derive(Debug, Clone, Serialize, Deserialize)]
401pub struct SimulationResult {
402    /// Time points
403    pub time: Vec<f64>,
404    /// Species concentrations over time
405    pub concentrations: HashMap<String, Vec<f64>>,
406    /// Reaction fluxes (optional)
407    pub fluxes: Option<HashMap<String, Vec<f64>>>,
408}
409
410/// COPASI-style simulator
411pub struct CopasiSimulation {
412    model: SbmlModel,
413    method: SimulationMethod,
414    /// Current state (concentrations)
415    state: Array1<f64>,
416    /// Current time
417    t: Time,
418    /// Time step
419    dt: Time,
420    /// RNG for stochastic simulations
421    rng_seed: u64,
422}
423
424impl CopasiSimulation {
425    /// Create new simulation
426    pub fn new(model: SbmlModel) -> Self {
427        let n = model.species.len();
428        let mut state = Array1::zeros(n);
429
430        // Initialize from model
431        for (i, species) in model.species.iter().enumerate() {
432            state[i] = species.initial_concentration.unwrap_or(0.0);
433        }
434
435        Self {
436            model,
437            method: SimulationMethod::Deterministic,
438            state,
439            t: 0.0,
440            dt: 0.01,
441            rng_seed: 42,
442        }
443    }
444
445    /// Set simulation method
446    pub fn set_method(&mut self, method: SimulationMethod) {
447        self.method = method;
448    }
449
450    /// Get current concentrations
451    pub fn get_concentrations(&self) -> HashMap<String, f64> {
452        self.model.species.iter()
453            .enumerate()
454            .map(|(i, s)| (s.id.clone(), self.state[i]))
455            .collect()
456    }
457
458    /// Run time course simulation
459    pub fn run(&mut self, duration: f64, n_points: usize) -> SimulationResult {
460        let dt = duration / n_points as f64;
461        let mut time = Vec::with_capacity(n_points + 1);
462        let mut concentrations: HashMap<String, Vec<f64>> = self.model.species.iter()
463            .map(|s| (s.id.clone(), Vec::with_capacity(n_points + 1)))
464            .collect();
465
466        // Record initial state
467        time.push(self.t);
468        for (i, species) in self.model.species.iter().enumerate() {
469            concentrations.get_mut(&species.id).unwrap().push(self.state[i]);
470        }
471
472        // Run simulation
473        for _ in 0..n_points {
474            self.step(dt);
475            time.push(self.t);
476            for (i, species) in self.model.species.iter().enumerate() {
477                concentrations.get_mut(&species.id).unwrap().push(self.state[i]);
478            }
479        }
480
481        SimulationResult {
482            time,
483            concentrations,
484            fluxes: None,
485        }
486    }
487
488    /// Single integration step
489    fn step(&mut self, dt: f64) {
490        match self.method {
491            SimulationMethod::Deterministic => self.step_deterministic(dt),
492            SimulationMethod::Stochastic => self.step_stochastic(),
493            SimulationMethod::TauLeaping => self.step_tau_leap(dt),
494            SimulationMethod::Hybrid => self.step_hybrid(dt),
495        }
496        self.t += dt;
497    }
498
499    /// Deterministic step (Euler, simple)
500    fn step_deterministic(&mut self, dt: f64) {
501        let rates = self.compute_rates();
502        let stoich = self.model.stoichiometry_matrix();
503
504        // dS/dt = N * v
505        let dstate = stoich.dot(&rates);
506        self.state = &self.state + &(&dstate * dt);
507
508        // Clamp to non-negative
509        for x in self.state.iter_mut() {
510            if *x < 0.0 {
511                *x = 0.0;
512            }
513        }
514    }
515
516    /// Stochastic step (Gillespie SSA)
517    fn step_stochastic(&mut self) {
518        // Simplified Gillespie - full implementation would use propensities
519        // and exponential waiting times
520        let rates = self.compute_rates();
521        let total_rate: f64 = rates.iter().sum();
522
523        if total_rate > 0.0 {
524            // Simple random selection (not proper SSA)
525            let dt = 1.0 / total_rate;
526            self.step_deterministic(dt);
527        }
528    }
529
530    /// Tau-leaping step
531    fn step_tau_leap(&mut self, tau: f64) {
532        // Simplified tau-leaping
533        self.step_deterministic(tau);
534    }
535
536    /// Hybrid step
537    fn step_hybrid(&mut self, dt: f64) {
538        // For now, just use deterministic
539        self.step_deterministic(dt);
540    }
541
542    /// Compute reaction rates
543    fn compute_rates(&self) -> Array1<f64> {
544        let n = self.model.reactions.len();
545        let mut rates = Array1::zeros(n);
546
547        for (j, reaction) in self.model.reactions.iter().enumerate() {
548            rates[j] = self.compute_reaction_rate(reaction);
549        }
550
551        rates
552    }
553
554    /// Compute rate for a single reaction
555    fn compute_reaction_rate(&self, reaction: &Reaction) -> f64 {
556        match &reaction.kinetic_law {
557            KineticLaw::MassAction { rate_constant } => {
558                let k = self.get_value(rate_constant);
559                let mut rate = k;
560                for sr in &reaction.reactants {
561                    let conc = self.get_species_concentration(&sr.species);
562                    rate *= conc.powf(sr.stoichiometry);
563                }
564                rate
565            }
566            KineticLaw::MichaelisMenten { vmax, km, substrate } => {
567                let vmax_val = self.get_value(vmax);
568                let km_val = self.get_value(km);
569                let s = self.get_species_concentration(substrate);
570                vmax_val * s / (km_val + s)
571            }
572            KineticLaw::Hill { vmax, k, substrate, n } => {
573                let vmax_val = self.get_value(vmax);
574                let k_val = self.get_value(k);
575                let s = self.get_species_concentration(substrate);
576                let s_n = s.powf(*n);
577                let k_n = k_val.powf(*n);
578                vmax_val * s_n / (k_n + s_n)
579            }
580            _ => 0.0,
581        }
582    }
583
584    /// Get parameter or species value
585    fn get_value(&self, id: &str) -> f64 {
586        // Try parameters first
587        if let Some(p) = self.model.get_parameter(id) {
588            return p.value;
589        }
590        // Then try species
591        self.get_species_concentration(id)
592    }
593
594    /// Get species concentration
595    fn get_species_concentration(&self, id: &str) -> f64 {
596        for (i, s) in self.model.species.iter().enumerate() {
597            if s.id == id {
598                return self.state[i];
599            }
600        }
601        0.0
602    }
603
604    /// Find steady state
605    pub fn steady_state(&mut self) -> Result<HashMap<String, f64>> {
606        // Simple iteration until convergence
607        let max_iter = 10000;
608        let tol = 1e-10;
609
610        for _ in 0..max_iter {
611            let old_state = self.state.clone();
612            self.step_deterministic(0.1);
613
614            let diff: f64 = (&self.state - &old_state)
615                .iter()
616                .map(|x| x.abs())
617                .sum();
618
619            if diff < tol {
620                return Ok(self.get_concentrations());
621            }
622        }
623
624        Err(OldiesError::NumericalError("Steady state not reached".into()))
625    }
626}
627
628// =============================================================================
629// STANDARD MODELS
630// =============================================================================
631
632pub mod models {
633    use super::*;
634
635    /// Michaelis-Menten enzyme kinetics
636    pub fn michaelis_menten() -> SbmlModel {
637        let mut model = SbmlModel::new("MichaelisMenten");
638
639        model.add_compartment(Compartment::new("cell", 1.0));
640        model.add_species(Species::new("S", "cell", 10.0));   // Substrate
641        model.add_species(Species::new("E", "cell", 1.0));    // Enzyme
642        model.add_species(Species::new("ES", "cell", 0.0));   // Complex
643        model.add_species(Species::new("P", "cell", 0.0));    // Product
644
645        model.add_parameter(Parameter::new("k1", 0.1));   // S + E -> ES
646        model.add_parameter(Parameter::new("k_1", 0.05)); // ES -> S + E
647        model.add_parameter(Parameter::new("k2", 0.1));   // ES -> E + P
648
649        // S + E -> ES
650        let mut r1 = Reaction::simple("binding", "S", "ES", "k1");
651        r1.reactants.push(SpeciesReference::new("E", 1.0));
652        model.add_reaction(r1);
653
654        // ES -> S + E
655        let mut r2 = Reaction::simple("unbinding", "ES", "S", "k_1");
656        r2.products.push(SpeciesReference::new("E", 1.0));
657        model.add_reaction(r2);
658
659        // ES -> E + P
660        let mut r3 = Reaction::simple("catalysis", "ES", "P", "k2");
661        r3.products.push(SpeciesReference::new("E", 1.0));
662        model.add_reaction(r3);
663
664        model
665    }
666
667    /// Lotka-Volterra predator-prey
668    pub fn lotka_volterra() -> SbmlModel {
669        let mut model = SbmlModel::new("LotkaVolterra");
670
671        model.add_compartment(Compartment::new("env", 1.0));
672        model.add_species(Species::new("prey", "env", 10.0));
673        model.add_species(Species::new("predator", "env", 5.0));
674
675        model.add_parameter(Parameter::new("alpha", 1.1));   // Prey growth
676        model.add_parameter(Parameter::new("beta", 0.4));    // Predation rate
677        model.add_parameter(Parameter::new("gamma", 0.4));   // Predator death
678        model.add_parameter(Parameter::new("delta", 0.1));   // Predator growth
679
680        model
681    }
682
683    /// Repressilator (synthetic gene network)
684    pub fn repressilator() -> SbmlModel {
685        let mut model = SbmlModel::new("Repressilator");
686
687        model.add_compartment(Compartment::new("cell", 1.0));
688
689        // mRNAs
690        model.add_species(Species::new("lacI", "cell", 0.0));
691        model.add_species(Species::new("tetR", "cell", 0.0));
692        model.add_species(Species::new("cI", "cell", 0.0));
693
694        // Proteins
695        model.add_species(Species::new("LacI", "cell", 0.0));
696        model.add_species(Species::new("TetR", "cell", 0.0));
697        model.add_species(Species::new("CI", "cell", 0.0));
698
699        model.add_parameter(Parameter::new("alpha", 216.0));
700        model.add_parameter(Parameter::new("alpha0", 0.216));
701        model.add_parameter(Parameter::new("beta", 5.0));
702        model.add_parameter(Parameter::new("n", 2.0));
703
704        model
705    }
706}
707
708// =============================================================================
709// TESTS
710// =============================================================================
711
712#[cfg(test)]
713mod tests {
714    use super::*;
715
716    #[test]
717    fn test_create_model() {
718        let model = models::michaelis_menten();
719        assert_eq!(model.species.len(), 4);
720        assert_eq!(model.reactions.len(), 3);
721    }
722
723    #[test]
724    fn test_stoichiometry_matrix() {
725        let model = models::michaelis_menten();
726        let stoich = model.stoichiometry_matrix();
727        assert_eq!(stoich.nrows(), 4);  // 4 species
728        assert_eq!(stoich.ncols(), 3);  // 3 reactions
729    }
730
731    #[test]
732    fn test_simulation() {
733        let model = models::michaelis_menten();
734        let mut sim = CopasiSimulation::new(model);
735        let result = sim.run(10.0, 100);
736
737        assert_eq!(result.time.len(), 101);
738        assert!(result.concentrations.contains_key("S"));
739        assert!(result.concentrations.contains_key("P"));
740    }
741
742    #[test]
743    fn test_mass_action_rate() {
744        let mut model = SbmlModel::new("test");
745        model.add_compartment(Compartment::new("c", 1.0));
746        model.add_species(Species::new("A", "c", 2.0));
747        model.add_species(Species::new("B", "c", 0.0));
748        model.add_parameter(Parameter::new("k", 0.1));
749        model.add_reaction(Reaction::simple("r1", "A", "B", "k"));
750
751        let sim = CopasiSimulation::new(model);
752        let conc = sim.get_concentrations();
753        assert_eq!(conc["A"], 2.0);
754    }
755}
oldies_copasi/lib.rs

oldies_copasi/
lib.rs
