means.approximation.mea package

Submodules

Gamma moment closure

This part of the package provides the original the Gamma closure.

class means.approximation.mea.closure_gamma.GammaClosure(max_order, multivariate=True)[source]

Bases: means.approximation.mea.closure_scalar.ClosureBase

EXPERIMENTAL

A class providing gamma closure to MomentExpansionApproximation. Expression for higher order (max_order + 1) central moments are computed from expressions of higher order raw moments. As a result, any higher order moments will be replaced by a symbolic expression depending on mean and variance only.

Parameters:
  • max_order (int) – the maximal order of moments to be modelled.
  • type – 0 for univariate (ignore covariances), 1 and 2 for

the two types of multivariate gamma distributions. :type type: int :return:

Log-normal moment closure

This part of the package provides the original the Log-normal closure.

class means.approximation.mea.closure_log_normal.LogNormalClosure(max_order, multivariate=True)[source]

Bases: means.approximation.mea.closure_scalar.ClosureBase

A class providing log-normal closure to MomentExpansionApproximation. Expression for higher order (max_order + 1) central moments are computed from expressions of higher order raw moments. As a result, any higher order moments will be replaced by a symbolic expression depending on mean and variance only.

Parameters:
  • max_order (int) – the maximal order of moments to be modelled.
  • multivariate – whether to consider covariances
Returns:

Normal moment closure

This part of the package provides the original the Normal (Gaussian) closure.

class means.approximation.mea.closure_normal.NormalClosure(max_order, multivariate=True)[source]

Bases: means.approximation.mea.closure_scalar.ClosureBase

A class providing normal closure to MomentExpansionApproximation. Expression for higher order (max_order + 1) central moments are directly computed using Isserlis’ Theorem. As a result, any higher order moments will be replaced by a symbolic expression depending on mean and variance only.

Parameters:
  • max_order (int) – the maximal order of moments to be modelled.
  • multivariate – whether to consider covariances
Returns:

Scalar moment closure

This part of the package provides the original (and default) closure ScalarClosure as well as the base class for all closers.

class means.approximation.mea.closure_scalar.ClosureBase(max_order, multivariate=True)[source]

Bases: object

A virtual class for closure methods. An implementation of _compute_raw_moments() must be provided in subclasses.

Parameters:
  • max_order (int) – the maximal order of moments to be modelled.
  • multivariate – whether to consider covariances
Returns:
close(mfk, central_from_raw_exprs, n_counter, k_counter)[source]

In MFK, replaces symbol for high order (order == max_order+1) by parametric expressions. That is expressions depending on lower order moments such as means, variances, covariances and so on.

Parameters:
  • mfk – the right hand side equations containing symbols for high order central moments
  • central_from_raw_exprs – expressions of central moments in terms of raw moments
  • n_counter (list[Moment]) – a list of Moments representing central moments
  • k_counter (list[Moment]) – a list of Moments representing raw moments
Returns:

the modified MFK
Return type:

sympy.Matrix
is_multivariate
max_order
class means.approximation.mea.closure_scalar.ScalarClosure(max_order, value=0)[source]

Bases: means.approximation.mea.closure_scalar.ClosureBase

A class providing scalar closure to MomentExpansionApproximation. Expression for higher order (max_order + 1) central moments are set to a scalar. Typically, higher order central moments are replaced by zero.

Parameters:
  • max_order (int) – the maximal order of moments to be modelled.
  • value – a scalar value for higher order moments
value
means.approximation.mea.dmu_over_dt.generate_dmu_over_dt(species, propensity, n_counter, stoichiometry_matrix)[source]

Calculate \(\frac{d\mu_i}{dt}\) in eq. 6 (see Ale et al. 2013).

\[\frac{d\mu_i}{dt} = S \begin{bmatrix} \sum_{l} \sum_{n_1=0}^{\infty} ... \sum_{n_d=0}^{\infty} \frac{1}{\mathbf{n!}} \frac{\partial^n \mathbf{n}a_l(\mathbf{x})}{\partial \mathbf{x^n}} |_{x=\mu} \mathbf{M_{x^n}} \end{bmatrix}\]
Parameters:
  • species (list[sympy.Symbol]) – the name of the species/variables (typically [‘y_0’, ‘y_1’, ..., ‘y_n’])
  • propensity – the reactions describes by the model
  • n_counter (list[Moment]) – a list of Moments representing central moments
  • stoichiometry_matrix (sympy.Matrix) – the stoichiometry matrix
Returns:

a matrix in which each row corresponds to a reaction, and each column to an element of counter.
means.approximation.mea.eq_central_moments.eq_central_moments(n_counter, k_counter, dmu_over_dt, species, propensities, stoichiometry_matrix, max_order)[source]

Function used to calculate the terms required for use in equations giving the time dependence of central moments.

The function returns the list Containing the sum of the following terms in in equation 9, for each of the \([n_1, ..., n_d]\) combinations in eq. 9 where ... is ... # FIXME

\[\mathbf{ {n \choose k} } (-1)^{ \mathbf{n-k} } [ \alpha \frac{d\beta}{dt} + \beta \frac{d\alpha}{dt} ]\]
Parameters:
  • n_counter (list[Moment]) – a list of Moments representing central moments
  • k_counter (list[Moment]) – a list of Moments representing raw moments
  • dmu_over_dt – du/dt in paper
  • species – species matrix: y_0, y_1,..., y_d
  • propensities – propensities matrix
  • stoichiometry_matrix – stoichiometry matrix
Returns:

central_moments matrix with (len(n_counter)-1) rows and one column per each \([n_1, ... n_d]\) combination
class means.approximation.mea.eq_mixed_moments.DBetaOverDtCalculator(propensities, n_counter, stoichoimetry_matrix, species)[source]

Bases: object

A class providing a efficient way to recursively calculate \(\frac{d\beta}{dt}\) (eq. 11 in [Ale2013]). A class was used here merely for optimisation reasons.

[Ale2013]
  1. Ale, P. Kirk, and M. P. H. Stumpf, “A general moment expansion method for stochastic kinetic models,” The Journal of Chemical Physics, vol. 138, no. 17, p. 174101, 2013.
Parameters:
  • propensities – the rates/propensities of the reactions
  • n_counter (list[Moment]) – a list of Moments representing central moments
  • stoichoimetry_matrix – The stoichiometry matrix. Explicitly provided by the model
  • species – the names of the variables/species
get(k_vec, e_counter)[source]

Provides the terms needed for equation 11 (see Ale et al. 2013). This gives the expressions for \(\frac{d\beta}{dt}\) in equation 9, these are the time dependencies of the mixed moments

Parameters:
  • k_vec – \(k\) in eq. 11
  • e_counter – \(e\) in eq. 11
Returns:

\(\frac{d\beta}{dt}\)

MEA helper functions.

This part of the package provides a few small utility functions for the rest mea.

means.approximation.mea.mea_helpers.derive_expr_from_counter_entry(expression, species, counter_entry)[source]

Derives an given expression with respect to arbitrary species and orders. This is used to compute \(\frac{\partial^n \mathbf{n}a_l(\mathbf{x})}{\partial \mathbf{x^n}}\) in eq. 6

Parameters:
  • expression (Expr) – the expression to be derived
  • species (list[Symbol]) – the name of the variables (typically {y_0, y_1, ..., y_n})
  • counter_entry – an entry of counter. That is a tuple of integers of length equal to the number of variables.

For example, (0,2,1) means we derive with respect to the third variable (first order) and to the second variable (second order)

Returns:the derived expression
means.approximation.mea.mea_helpers.get_one_over_n_factorial(*args)[source]

Calculates the \(\frac{1}{\mathbf{n!}}\) of eq. 6 (see Ale et al. 2013). That is the invert of a product of factorials. :param counter_entry: an entry of counter. That is an array of integers of length equal to the number of variables. For instance, counter_entry could be [1,0,1] for three variables. :return: a scalar as a sympy expression

means.approximation.mea.mea_helpers.make_k_chose_e(e_vec, k_vec)[source]

Computes the product \({\mathbf{n} \choose \mathbf{k}}\)

Parameters:
  • e_vec (numpy.array) – the vector e
  • k_vec (numpy.array) – the vector k
Returns:

a scalar
class means.approximation.mea.moment_expansion_approximation.MomentExpansionApproximation(model, max_order, closure='scalar', *closure_args, **closure_kwargs)[source]

Bases: means.approximation.approximation_baseclass.ApproximationBaseClass

A class to perform moment expansion approximation as described in [Ale2013] up to a given order of moment. In addition, it allows to close the Taylor expansion by using parametric values for last order central moments.

Parameters:
  • model (Model) – The model to be approximated
  • max_order – the highest order of central moments in the resulting ODEs
  • closure (string) –

    a string describing the type of closure to use. Currently, the supported closures are:

    ‘scalar’
    higher order central moments are set to zero. See ScalarClosure.
    ‘normal’
    uses normal distribution to compute last order central moments. See NormalClosure.
    ‘log-normal’
    uses log-normal distribution. See LogNormalClosure.
    ‘gamma’
    EXPERIMENTAL, uses gamma distribution. See GammaClosure.
  • closure_args – arguments to be passed to the closure
  • closure_kwargs – keyword arguments to be passed to the closure
closure
run()[source]

Overrides the default run() method. Performs the complete analysis on the model specified during initialisation.

Returns:an ODE problem which can be further used in inference and simulation.
Return type:ODEProblem
means.approximation.mea.moment_expansion_approximation.mea_approximation(model, max_order, closure='scalar', *closure_args, **closure_kwargs)[source]

A wrapper around MomentExpansionApproximation. It performs moment expansion approximation (MEA) up to a given order of moment. See MomentExpansionApproximation for details about the options.

Returns:an ODE problem which can be further used in inference and simulation.
Return type:ODEProblem
means.approximation.mea.raw_to_central.raw_to_central(n_counter, species, k_counter)[source]

Expresses central moments in terms of raw moments (and other central moments). Based on equation 8 in the paper:

\[\mathbf{M_{x^n}} = \sum_{k_1=0}^{n_1} ... \sum_{k_d=0}^{n_d} \mathbf{{n \choose k}} (-1)^{\mathbf{n-k}} \mu^{\mathbf{n-k}} \langle \mathbf{x^k} \rangle\]

The term \(\mu^{\mathbf{n-k}}\), so called alpha term is expressed with respect to species values that are equivalent to \(\mu_i\) in the paper.

The last term, the beta term, \(\langle \mathbf{x^n} \rangle\) is simply obtained from k_counter as it contains the symbols for raw moments.

Parameters:
  • n_counter (list[Moment]) – a list of Moments representing central moments
  • species – the symbols for species means
  • k_counter (list[Moment]) – a list of Moments representing raw moments
Returns:

a vector of central moments expressed in terms of raw moment

Module contents

Moment Expansion Approximation

This part of the package implements Moment Expansion Approximation as described in [Ale2013]. In addition to the standard implementation, it allows to use different distribution (such as normal, log-normal and gamma) to close the moment expansion. The function mea_approximation() should provide all the necessary options.

Example:
>>> from means import mea_approximation
>>> from means.examples.sample_models import MODEL_P53
>>> ode_problem = mea_approximation(MODEL_P53,max_order=2)
>>> # equivalent to
>>> # ode_problem = mea_approximation(MODEL_P53, max_order=2, closure="scalar", value=0)
>>> print ode_problem

The result is an means.core.problems.ODEProblem. Typically, it would be further used to perform simulations (see simulation) and inference (see inference).

[Ale2013]
  1. Ale, P. Kirk, and M. P. H. Stumpf, “A general moment expansion method for stochastic kinetic models,” The Journal of Chemical Physics, vol. 138, no. 17, p. 174101, 2013.