crn 0.1.2

Stochastic and deterministic simulation of chemical reaction networks.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121

use rand::Rng;

use crate::{state::State, Crn};

use std::fmt::Debug;

use thiserror::Error;

/// A simulation can fail because no more reactions are possible, or because of numerical instability.
#[derive(Error, Debug)]
pub enum Error {
    /// No reactions are possible from the current state.
    #[error("CRN has reached terminal state")]
    TerminalState,
    /// The simulation has become numerically unstable.
    #[error("Insufficient precision for accurate simulation")]
    InsufficientPrecision,
}

/// A stochastic CRN. This is simulated using the Gillespie algorithm. Stochastic CRNs are essentially a type of continuous-time Markov chain.
pub type StoCrn = Crn<i32>;

impl StoCrn {
    /// Simulate one reaction. Uses `rates` to avoid repeated allocations.
    fn step(&mut self, rates: &mut [f64]) -> Result<(), Error> {
        let mut rate = 0.0;

        self.rxns.iter().enumerate().for_each(|(idx, rxn)| {
            let cur_rate = self.state.rate(rxn);
            rates[idx] = cur_rate;
            rate += cur_rate;
        });

        if rate == 0.0 {
            return Err(Error::TerminalState);
        }

        let mut rng = rand::thread_rng();
        // the random number is in (0, 1], so the ln is negative or zero and this is really an addition
        self.state.time -= (1.0 - rng.gen::<f64>()).ln() / rate;
        let j = rng.gen::<f64>() * rate;
        let mut sum = 0.0;

        for (idx, cur_rate) in rates.iter().enumerate() {
            sum += cur_rate;
            if j < sum {
                self.state.apply(&self.rxns[idx]);
                return Ok(());
            }
        }
        Err(Error::InsufficientPrecision)
    }

    /// Simulate a number of reactions.
    pub fn steps(&mut self, steps: usize) -> Result<(), Error> {
        let mut rates = vec![0.0; self.rxns.len()];
        for _ in 0..steps {
            self.step(&mut rates)?;
        }
        Ok(())
    }

    /// Simulates for a given amount of time. Returns a collection of individual species' history.
    pub fn simulate_history(&mut self, t: f64) -> Result<Vec<State<f64>>, Error> {
        let mut result = Vec::new();

        let mut rates = vec![0.0; self.rxns.len()];
        while self.state.time < t {
            if self.step(&mut rates).is_err() {
                break;
            }
            let species = self.state.species.iter().map(|x| *x as f64).collect();
            result.push(State {
                species,
                time: self.state.time,
            });
        }
        Ok(result)
    }
}

#[cfg(test)]
mod tests {
    use crate::StoCrn;

    #[test]
    fn test() {
        const N: i32 = 100;
        let mut crn = StoCrn::parse(&format!("A = {N}; A -> ;")).unwrap();
        let mut rates = vec![0.0; crn.rxns.len()];
        for i in (0..N).rev() {
            crn.step(&mut rates).unwrap();
            assert_eq!(crn.state.species[0], i);
        }
    }

    #[test]
    fn test2() {
        let mut crn = StoCrn::parse("A = 1; B = 1; A + B -> C; C -> A + B;").unwrap();
        let mut rates = vec![0.0; crn.rxns.len()];
        crn.step(&mut rates).unwrap();
        assert_eq!(crn.state.species[0], 0);
        assert_eq!(crn.state.species[1], 0);
        assert_eq!(crn.state.species[2], 1);
        crn.step(&mut rates).unwrap();
        assert_eq!(crn.state.species[0], 1);
        assert_eq!(crn.state.species[1], 1);
        assert_eq!(crn.state.species[2], 0);
    }

    #[test]
    fn test3() {
        let mut crn = StoCrn::parse("A = 16; 2A -> B; 2B -> C; 2C -> D; 2D -> E;").unwrap();
        crn.steps(15).unwrap();
        assert_eq!(crn.state.species[0], 0);
        assert_eq!(crn.state.species[1], 0);
        assert_eq!(crn.state.species[2], 0);
        assert_eq!(crn.state.species[3], 0);
        assert_eq!(crn.state.species[4], 1);
    }
}