diff --git a/arrow/arrow.py b/arrow/arrow.py index a7d3183..8ecae92 100644 --- a/arrow/arrow.py +++ b/arrow/arrow.py @@ -80,12 +80,13 @@ class StochasticSystem(object): (and zero everywhere else). ''' - def __init__(self, stoichiometry, rates, random_seed=0): + def __init__(self, stoichiometry, rates, forms=None, random_seed=0): ''' This invokes the Obsidian C module (see obsidianmodule.c) with the - stoichiometry, reaction rates and a variety of derived values. Once constructed, - this can be invoked by calling `evolve` with a duration and initial state, since - the stoichiometry will be shared among all calls. + stoichiometry, reaction rates, the form of each reaction (optional) + and a variety of derived values. Once constructed, this can be invoked + by calling `evolve` with a duration and initial state, since the + stoichiometry will be shared among all calls. There are four derived values, each of which is a list of variable length lists. In order to pass this into C, these nested lists are flattened and two @@ -101,7 +102,14 @@ def __init__(self, stoichiometry, rates, random_seed=0): ''' self.stoichiometry = stoichiometry - self.rates = rates + if rates.ndim == 1: + rates = [[rate] for rate in rates.astype(np.float64)] + self.rates_flat, self.rates_lengths, self.rates_indexes = flat_indexes(rates) + + if forms is not None: + self.forms = forms + else: + self.forms = np.zeros(stoichiometry.shape[0], dtype=np.int64) self.random_seed = random_seed @@ -117,7 +125,10 @@ def __init__(self, stoichiometry, rates, random_seed=0): self.obsidian = obsidian.obsidian( self.random_seed, self.stoichiometry, - self.rates, + self.rates_flat, + self.rates_lengths, + self.rates_indexes, + self.forms, self.reactants_lengths, self.reactants_indexes, self.reactants_flat, @@ -131,7 +142,7 @@ def __init__(self, stoichiometry, rates, random_seed=0): def evolve(self, duration, state): steps, time, events, outcome = self.obsidian.evolve(duration, state) - occurrences = np.zeros(len(self.rates)) + occurrences = np.zeros(self.stoichiometry.shape[0]) for event in events: occurrences[event] += 1 diff --git a/arrow/obsidian.c b/arrow/obsidian.c index ab77a33..0270bfd 100644 --- a/arrow/obsidian.c +++ b/arrow/obsidian.c @@ -24,6 +24,12 @@ choose(long n, long k) { return combinations; } +// Compute the fractional saturation +double +fraction_saturation(long s, double k) { + return ((double) s / (s + k)); +} + // Perform the Gillespie algorithm with the given stoichiometry, reaction rates and // initial state for the provided duration. @@ -34,11 +40,14 @@ choose(long n, long k) { // contains all the information about the reactions this system performs. Each // row is a reaction, and each column is a substrate, so every entry is the // change in counts of each substrate when the given reaction is performed. -// * rates: an array of length `reactions_count` that encodes the base rate for -// each reaction. The actual propensity is further dependent on the counts for -// each reactant in the reaction. +// * rates_flat: an array of length `reactions_count` (or greater if the reaction +// is not of the default form 0) that encodes the base rate for each reaction. +// The actual propensity is further dependent on the counts for each reactant +// in the reaction. +// * forms: an array of length `reactions_count` that indicates the type of +// propensity computation that is associated with each reaction. // * duration: How long to run the simulation for. -// * state: An array of length `substrates_count` that contains the inital count +// * state: An array of length `substrates_count` that contains the initial count // for each substrate. The outcome of the algorithm will involve an array of // the same size that signifies the counts of each substrate after all the // reactions are performed. @@ -83,7 +92,10 @@ evolve(MTState *random_state, int reactions_count, int substrates_count, long *stoichiometry, - double *rates, + double *rates_flat, + long *rates_lengths, + long *rates_indexes, + long *forms, long *reactants_lengths, long *reactants_indexes, @@ -158,18 +170,47 @@ evolve(MTState *random_state, // First update all propensities that were affected by the previous event for (up = 0; up < update_length; up++) { - // Find which reaction to update and initialize the propensity for that reaction - // with the rate for that reaction + // Find which reaction to update reaction = update[up]; - propensities[reaction] = rates[reaction]; - - // Go through each reactant and calculate its contribution to the propensity - // based on the counts of the corresponding substrate, which is then multiplied - // by the reaction's original rate and the contributions from other reactants - for (reactant = 0; reactant < reactants_lengths[reaction]; reactant++) { - index = reactants_indexes[reaction] + reactant; - count = outcome[reactants[index]]; - propensities[reaction] *= choose(count, reactions[index]); + + // Compute propensity depending on the reaction's form + switch(forms[reaction]) { + + case 0: // Standard Gillespie + // Initialize the propensity for reaction with its rate constant + propensities[reaction] = rates_flat[reaction]; + + // Go through each reactant and calculate its contribution to the + // propensity based on the counts of the corresponding substrate, + // which is then multiplied by the reaction's original rate and the + // contributions from other reactants + for (reactant = 0; reactant < reactants_lengths[reaction]; reactant++) { + index = reactants_indexes[reaction] + reactant; + count = outcome[reactants[index]]; + propensities[reaction] *= choose(count, reactions[index]); + } + + break; + + case 1: // Michaelis-Menten + // Initialize the propensity for reaction with its maximal propensity + index = rates_indexes[reaction]; + propensities[reaction] = rates_flat[index]; + + // Go through each reactant and calculate its contribution to the + // propensity based on the counts of the corresponding substrates + // relative to their Michaelis-Menten constant. + for (reactant = 0; reactant < reactants_lengths[reaction]; reactant++) { + index = reactants_indexes[reaction] + reactant; + count = outcome[reactants[index]]; + propensities[reaction] *= fraction_saturation(count, + rates_flat[rates_indexes[reaction] + reactant + 1]); + } + + break; + + default: + printf("arrow.obsidian.evolve - unexpected form: %ld", forms[reaction]); } } diff --git a/arrow/obsidian.h b/arrow/obsidian.h index ad70f6a..ab06354 100644 --- a/arrow/obsidian.h +++ b/arrow/obsidian.h @@ -19,7 +19,10 @@ evolve(MTState *random_state, int reactions_count, int substrates_count, long *stoichiometry, - double *rates, + double *rates_flat, + long *rates_lengths, + long *rates_indexes, + long *forms, long *reactants_lengths, long *reactants_indexes, diff --git a/arrow/obsidianmodule.c b/arrow/obsidianmodule.c index 4e4ee1d..2f310fa 100644 --- a/arrow/obsidianmodule.c +++ b/arrow/obsidianmodule.c @@ -39,7 +39,10 @@ typedef struct { int reactions_count; int substrates_count; long *stoichiometry; - double *rates; + double *rates_flat; + long *rates_lengths; + long *rates_indexes; + long *forms; long *reactants_lengths; long *reactants_indexes; @@ -66,7 +69,10 @@ newObsidianObject(int random_seed, int reactions_count, int substrates_count, long *stoichiometry, - double *rates, + double *rates_flat, + long *rates_lengths, + long *rates_indexes, + long *forms, long *reactants_lengths, long *reactants_indexes, @@ -103,7 +109,10 @@ newObsidianObject(int random_seed, self->reactions_count = reactions_count; self->substrates_count = substrates_count; self->stoichiometry = stoichiometry; - self->rates = rates; + self->rates_flat = rates_flat; + self->rates_lengths = rates_lengths; + self->rates_indexes = rates_indexes; + self->forms = forms; self->reactants_lengths = reactants_lengths; self->reactants_indexes = reactants_indexes; @@ -138,7 +147,7 @@ Obsidian_demo(ObsidianObject *self, PyObject *args) if (!PyArg_ParseTuple(args, ":demo")) return NULL; - print_array(self->rates, self->reactions_count); + print_array(self->rates_flat, self->reactions_count); Py_INCREF(Py_None); return Py_None; @@ -187,7 +196,10 @@ Obsidian_evolve(ObsidianObject *self, PyObject *args) self->reactions_count, self->substrates_count, self->stoichiometry, - self->rates, + self->rates_flat, + self->rates_lengths, + self->rates_indexes, + self->forms, self->reactants_lengths, self->reactants_indexes, self->reactants, @@ -364,7 +376,10 @@ _invoke_obsidian(PyObject *self, PyObject *args) { int random_seed; PyObject *stoichiometry_obj, - *rates_obj, + *rates_flat_obj, + *rates_lengths_obj, + *rates_indexes_obj, + *forms_obj, *reactants_lengths_obj, *reactants_indexes_obj, @@ -380,10 +395,13 @@ _invoke_obsidian(PyObject *self, PyObject *args) { *substrates_obj; if (!PyArg_ParseTuple(args, - "iOOOOOOOOOOOO", + "iOOOOOOOOOOOOOOO", &random_seed, &stoichiometry_obj, - &rates_obj, + &rates_flat_obj, + &rates_lengths_obj, + &rates_indexes_obj, + &forms_obj, &reactants_lengths_obj, &reactants_indexes_obj, @@ -406,8 +424,16 @@ _invoke_obsidian(PyObject *self, PyObject *args) { long *stoichiometry = (long *) PyArray_DATA(stoichiometry_array); // Import the rates for each reaction - PyObject *rates_array = array_for(rates_obj, NPY_DOUBLE); - double *rates = (double *) PyArray_DATA(rates_array); + PyObject *rates_flat_array = array_for(rates_flat_obj, NPY_DOUBLE); + double *rates_flat = (double *) PyArray_DATA(rates_flat_array); + PyObject *rates_lengths_array = array_for(rates_lengths_obj, NPY_INT64); + long *rates_lengths = (long *) PyArray_DATA(rates_lengths_array); + PyObject *rates_indexes_array = array_for(rates_indexes_obj, NPY_INT64); + long *rates_indexes = (long *) PyArray_DATA(rates_indexes_array); + + // Import the forms for each reaction + PyObject *forms_array = array_for(forms_obj, NPY_INT64); + long *forms = (long *) PyArray_DATA(forms_array); // Pull all the precalculated nested arrays PyObject *reactants_lengths_array = array_for(reactants_lengths_obj, NPY_INT64); @@ -438,7 +464,10 @@ _invoke_obsidian(PyObject *self, PyObject *args) { reactions_count, substrates_count, stoichiometry, - rates, + rates_flat, + rates_lengths, + rates_indexes, + forms, reactants_lengths, reactants_indexes, @@ -459,6 +488,10 @@ _invoke_obsidian(PyObject *self, PyObject *args) { } // Clean up all the PyObject * references + Py_XDECREF(rates_flat_array); + Py_XDECREF(rates_lengths_array); + Py_XDECREF(rates_indexes_array); + Py_XDECREF(reactants_lengths_array); Py_XDECREF(reactants_indexes_array); Py_XDECREF(reactants_array); diff --git a/arrow/test/test_arrow.py b/arrow/test/test_arrow.py index b2f9b83..6f1f130 100644 --- a/arrow/test/test_arrow.py +++ b/arrow/test/test_arrow.py @@ -173,6 +173,25 @@ def test_compare_runtime(): print('reference time elapsed: {}'.format(reference_end - reference_start)) print('obsidian time elapsed: {}'.format(obsidian_end - obsidian_start)) +def test_forms(): + stoichiometric_matrix = np.array([ + [-1, -1, 1, 0, 0, 0], + [0, 0, 0, -1, -1, 1]], np.int64) + + rates = np.array([ + [10, 1, 1], + [1, 1, 1]], np.float64) + forms = np.array([1, 1]) + + arrow = StochasticSystem(stoichiometric_matrix, rates, forms) + result = arrow.evolve(1., np.array([50, 50, 0, 20, 20, 0], np.int64)) + + print('steps: {}'.format(result['steps'])) + print('events: {}'.format(result['events'])) + print('occurrences: {}'.format(result['occurrences'])) + print('outcome: {}'.format(result['outcome'])) + + if __name__ == '__main__': parser = argparse.ArgumentParser() parser.add_argument('--plot', action='store_true')