Skip to content

gpac API

This is the API documentation for gpac. See the Github page for examples of usage and installation instructions.

gpac is a Python package for numerically simulating a general-purpose analog computer (GPAC), defined by Claude Shannon in 1941 as an abstract model of programmable analog computational devices such as the differential analyzer created by Vannevar Bush and Harold Locke Hazen in the 1920s. General descriptions of GPACs can be found here and here.

It also has support for a very common model governed by polynomial ODEs, that of continuous mass-action chemical reaction networks. And despite having nothing to do with GPAC or ODEs, it also can simulate discrete CRNs using the excellent rebop package; see the functions plot_gillespie, rebop_crn_counts, and rebop_sample_future_configurations.

GPACs are typically defined by a circuit with gates that can add, multiply, introduce constants, and integrate an input with respect to time. The most elegant way to specify a GPAC is by defining a set of ordinary differential equations (ODEs) corresponding to the output wires of integrator gates in the GPAC circuit.

So essentially, this package makes it easy to write down such ODEs and numerically integrate and plot them.

Although gpac has two submodules ode and crn, you can import all elements from both directly from gpac, e.g., from gpac import plot, plot_crn.

There are many examples in the jupyter notebook on Github.

API Reference

gpac.ode

This module has ODE-only functions, mainly integrate_odes and plot.

ValOde module-attribute 

ValOde = TypeVar('ValOde', 'sympy.Expr', float, int)

A type variable representing a value that can be a sympy expression, float, or int. This represents the type of the values in the odes dict passed to integrate_odes and plot and similar functions such as plot_crn. The "typical" case is a sympy expression such as y - x*y in odes = { x: y - x*y }, where x and y are sympy symbols. However, we may wish to have a constant derivative such as odes = { t: 1 }; this avoids the user having to wrap the constant 1 in a sympy expression via sympy.sympify(1).

Config module-attribute 

Config: TypeAlias = Mapping['sympy.Symbol', float]

Type alias for a configuration, such as the inits parameter of integrate_odes and plot representing the initial configuration of the system. It is a dict mapping each variable to its value.

Number module-attribute 

Number = TypeVar('Number', int, float)

Type variable representing a number, either an int or a float. This is required to avoid some type hint errors due to the fact that Mapping is invariant in its key type, so for the resets parameter of some functions, we declare it to be type Mapping[Number, Config] rather than Mapping[float, Config], which would cause a type checker error trying to declare a resets dict with int keys such as resets = {1: {a: 4.5}, 2: {b: 6.5}}.

default_figsize module-attribute 

default_figsize = (12, 3)

Default figure size for matplotlib plotting functions such as plot and plot_crn.

integrate_odes 

integrate_odes(
    odes: Mapping[Symbol, ValOde],
    inits: Config,
    t_eval: Iterable[float] | None = None,
    *,
    t_span: tuple[float, float] | None = None,
    dependent_symbols: Iterable[ValOde] = (),
    resets: Mapping[Number, Config] | None = None,
    method: str | OdeSolver = "RK45",
    dense_output: bool = False,
    events: Callable | Iterable[Callable] | None = None,
    vectorized: bool = False,
    args: tuple | None = None,
    **options
) -> OdeResult

Integrate the given ODEs using scipy's solve_ivp function in the package scipy.integrate, returning the same object returned by solve_ivp.

This is a convienence function that wraps the scipy function solve_ivp, allowing the user to specify the ODEs using sympy symbols and expressions (instead of a Python function on tuples of floats, which is what solve_ivp expects, but is more awkward to specify than using sympy expressions).

The object solution returned by solve_ivp has field solution.y which is a 2D numpy array, each row of which is the trajectory of a value in the ODEs. The order of the rows is the same as the iteration order of the keys in the odes dict.

Besides the parameters described below, all other parameters are simply passed along to solve_ivp in scipy.integrate. As with that function, the following are explicitly named parameters: method, dense_output, events, vectorized, args, and all other keyword arguments are passed in through **options; see the documentation for solve_ivp for a description of these parameters: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

rock-paper-scissors oscillator example
import sympy, gpac, numpy as np
a,b,c = sympy.symbols('a b c')
odes = {
    a: -a*b + c*a,
    b: -b*c + a*b,
    c: -c*a + b*c,
}
inits = {
    a: 10,
    b: 1,
    c: 1,
}
t_eval = np.linspace(0, 1, 5)
print(gpac.integrate_odes(odes, inits, t_eval=t_eval))

This outputs 

 message: 'The solver successfully reached the end of the integration interval.'
    nfev: 62
    njev: 0
     nlu: 0
     sol: None
  status: 0
 success: True
       t: array([0.  , 0.25, 0.5 , 0.75, 1.  ])
t_events: None
       y: array([[10.        ,  4.84701622,  0.58753815,  0.38765743,  3.07392998],
      [ 1.        ,  6.84903338,  9.63512628,  3.03634559,  0.38421121],
      [ 1.        ,  0.3039504 ,  1.77733557,  8.57599698,  8.54185881]])
y_events: None

All symbols are interpreted as functions of a single variable called "time", and the derivatives are with respective to time.

Although you cannot reference the time variable directly in the ODEs, this can be simulated by introducing a new variable t whose derivative is 1 and initial value is the initial time. For example, the following code implements a(t) = sin(t) (with time derivative a'(t) = cos(t)) and b(t) = -(t/2 - 1)^2 + 2 (with time derivative b'(t) = 1 - t/2). In this case we want t=0 initially, and any symbol that is not a key in inits is assumed to have initial value 0.

# trick for referencing time variable directly in ODEs
from sympy import sin, cos
from math import pi

a,b,t = sympy.symbols('a b t')
odes = {
    a: cos(t),
    b: 1 - t/2, # derivative of -(t/2 - 1)^2 + 2
    t: 1,
}
inits = { b: 1 }
t_eval = np.linspace(0, 2*pi, 5) # [0, pi/2, pi, 3*pi/2, 2*pi]
solution = gpac.integrate_odes(odes, inits, t_eval=t_eval)
print(f'a(pi/2) = {solution.y[0][1]:.2f}')
print(f'a(pi)   = {solution.y[0][2]:.2f}')
print(f'b(pi/2) = {solution.y[1][1]:.2f}')
print(f'b(pi)   = {solution.y[1][2]:.2f}')

which prints

a(pi/2) = 1.00
a(pi)   = 0.00
b(pi/2) = 1.95
b(pi)   = 1.67

Parameters:

Name Type Description Default
odes Mapping[Symbol, ValOde]

dict mapping sympy symbols to sympy expressions representing the ODEs. Alternatively, the keys can be strings, and the values can be strings that look like expressions, e.g., {'a': '-a*b + c*a'}. If a symbol is referenced in an expression but is not a key in odes, a ValueError is raised.

 required
inits Config

dict mapping sympy symbols to initial values of each symbol. Alternatively, the keys can be strings. Any symbols in the ODEs that are not keys in inits will be assumed to have initial value of 0. If a symbol appears as a key in inits but is not a key in odes, a ValueError is raised.

 required
t_eval Iterable[float] | None

iterable of times at which to evaluate the ODEs. At least one of t_eval or t_span must be specified.

 None
t_span tuple[float, float] | None

pair (start_time, end_time) for the integration. If not specified, first and last times in t_eval are used. (This is different from solve_ivp, which requires t_span to be specified.) At least one of t_eval or t_span must be specified.

 None
dependent_symbols Iterable[ValOde]

iterable of sympy expressions (or strings) representing symbols that are functions of the other symbols that are keys in odes. These values are added to the end of the 2D array field sol.y in the object sol returned by solve_ivp, in the order in which they appear in dependent_variables. For an example, see the example notebook. Note that this is a different type than the dependent_symbols parameter in plot, which is a dict mapping sympy Symbols to sympy expressions. In other words in this function, the dependent symbols have no "names", only positions, corresponding to the fact that solve_ivp does not have named symbols, only positions within the vector of solutions. There is no point in giving names to the dependent symbols here, since the names are lost after handing the ODEs to solve_ivp. However, in plot and similar functions such as plot_crn, the names are used to label the plot.

 ()
resets Mapping[Number, Config] | None

If specified, this is a dict mapping times to "configurations" (i.e., dict mapping symbols to values). The configurations are used to set the values of the symbols manually during the ODE integration at specific times. Any symbols not appearing as keys in resets are left at their current values. The OdeResult returned (the one returned by solve_ivp in scipy) will have two additional fields: reset_times and reset_indices, which are lists of the times and indices in sol.t corresponding to the times when the resets were applied. Raises a ValueError if any time lies outside the integration interval, or if resets is empty.

 None
method str | OdeSolver

See solve_ivp

 'RK45'
dense_output bool

See solve_ivp

 False
events Callable | Iterable[Callable] | None

See solve_ivp

 None
vectorized bool

See solve_ivp

 False
args tuple | None

See solve_ivp

 None
options

This is a catch-all for any additional keyword arguments that are passed to solve_ivp, for example you could pass rtol=1e-6 to set the relative tolerance to 1e-6:

plot(odes, inits, t_eval=t_eval, rtol=1e-6)

For solver-specific parameters, see solve_ivp

 {}

Returns:

Type Description
OdeResult

solution to the ODEs, same as object returned by solve_ivp

plot 

plot(
    odes: dict[Symbol, ValOde],
    inits: Config,
    t_eval: Iterable[float] | None = None,
    *,
    t_span: tuple[float, float] | None = None,
    resets: Mapping[Number, Config] | None = None,
    dependent_symbols: dict[Symbol, ValOde] | None = None,
    figsize: tuple[float, float] = default_figsize,
    symbols_to_plot: (
        Iterable[Symbol]
        | Iterable[Sequence[Symbol]]
        | str
        | Pattern
        | Iterable[Pattern]
        | None
    ) = None,
    legend: dict[Symbol, str] | None = None,
    latex_legend: bool = False,
    omit_legend: bool = False,
    show: bool = False,
    method: str | OdeSolver = "RK45",
    dense_output: bool = False,
    events: Callable | Iterable[Callable] | None = None,
    vectorized: bool = False,
    return_ode_result: bool = False,
    args: tuple | None = None,
    loc: str | tuple[float, float] = "best",
    warn_change_dpi: bool = False,
    **options
) -> OdeResult | None

Numerically integrate the given ODEs using the function integrate_odes, then plot the trajectories using matplotlib. (Assumes it is being run in a Jupyter notebook.)

Parameters:

Name Type Description Default
odes dict[Symbol, ValOde]

See integrate_odes.

 required
inits Config

See integrate_odes.

 required
t_eval Iterable[float] | None

See integrate_odes.

 None
t_span tuple[float, float] | None

See integrate_odes.

 None
dependent_symbols dict[Symbol, ValOde] | None

See integrate_odes.

 None
resets Mapping[Number, Config] | None

See integrate_odes.

 None
method str | OdeSolver

See solve_ivp

 'RK45'
dense_output bool

See solve_ivp

 False
events Callable | Iterable[Callable] | None

See solve_ivp

 None
vectorized bool

See solve_ivp

 False
args tuple | None

See solve_ivp

 None
figsize tuple[float, float]

pair (width, height) of the figure

 default_figsize
latex_legend bool

If True, surround each symbol name with dollar signs, unless it is already surrounded with them, so that the legend is interpreted as LaTeX. If this is True, then the symbol name must either start and end with $, or neither start nor end with $.

 False
symbols_to_plot Iterable[Symbol] | Iterable[Sequence[Symbol]] | str | Pattern | Iterable[Pattern] | None

symbols to plot; if not specified, then all symbols are plotted. If it is a 2D list (or other Iterable of Iterables of strings or symbols), then each group of symbols is plotted in a separate subplot. If a string or re.Pattern, then only symbols whose names match the string or pattern are plotted.

 None
legend dict[Symbol, str] | None

If specified, should be a dict mapping symbols (or strings) to strings. For each symbol that is plotted, the corresponding string is used as the label in the plot's legend instead of the original name of the symbol. This can be useful for example to include LaTeX, mapping a symbol with a name like 'xt' to a string like r'$x_t$'.

 None
omit_legend bool

If True, do not show the legend at all. Raises exception if true and legend is specified or latex_legend is also true.

 False
show bool

whether to call matplotlib.pyplot.show() after creating the plot; If False, this helps the user to call other functions such as matplotlib.pyplot.legend() or matplotlib.pyplot.grid() after calling this function, which will not work if matplotlib.pyplot.show() has already been called. However, if you want to display multiple plots from the same cell in a Jupyter notebook, you should either set this to True, or (in case you want to configure each plot by calling other matplotlib.pyplot functions, such as yscale), manually call matplotlib.pyplot.show() after each call to this function.

 False
dependent_symbols dict[Symbol, ValOde] | None

dict mapping symbols (or strings) to sympy expressions (or strings) representing variables that are functions of the other variables that are keys in odes. For an example, see the [example notebook] (https://github.com/UC-Davis-molecular-computing/gpac/blob/main/notebook.ipynb). Note that this is a different type than the dependent_symbols parameter in integrate_odes, which is an Iterable of sympy expressions.

 None
return_ode_result bool

if True, returns solution to the ODEs, same as object returned by solve_ivp in scipy.integrate. Otherwise (default) None is returned. The reason the solution is not automatically returned is that it pollutes the output of a Jupyter notebook, so this avoids needing to type something like _ = gpac.plot(...) to suppress the output. But if you want that solution object, you can set this to True.

 False
loc str | tuple[float, float]

location of the legend; see documentation for matplotlib.pyplot.legend

 'best'
warn_change_dpi bool

If True, print a warning if the dpi of the figure gets changed from its default.

 False
options

See integrate_odes.

 {}

Returns:

Type Description
OdeResult | None

Typically None, but if return_ode_result is True, returns the solution to the ODEs, same as object returned by solve_ivp.

display_odes 

display_odes(odes: dict[Symbol, ValOde]) -> None

Display the ODEs in a readable format in a Jupyter notebook.

Parameters:

Name Type Description Default
odes dict[Symbol, ValOde]

dict mapping sympy symbols (or strings) to sympy expressions (or strings or floats) representing the ODEs. Alternatively, the keys can be strings, and the values can be strings that look like expressions, e.g., {'a': '-a*b + c*a'}.

 required

gpac.crn

Module for expressing chemical reaction networks and deriving their ODEs. Ideas and much code taken from this repo.

For example, to specify the "approximate majority" chemical reaction network (see DOI: 10.1007/978-3-540-75142-7_5 or DOI: 10.1126/science.aal2052)

\[ \begin{align*} A+B &\to 2U \\ A+U &\to 2A \\ B+U &\to 2B \end{align*} \]

we can write

import gpac
import numpy as np

a, b, u = gpac.species('A B U')
rxns = [
    a+b >> 2*u,
    a+u >> 2*a,
    b+u >> 2*b,
]
inits = {a: 0.51, b: 0.49}
t_eval = np.linspace(0, 10, 100)
gpac.plot_crn(rxns, inits, t_eval)

which will plot the concentrations of A, B, and U over time. One can specify reversible reactions by using the | operator instead of >> (e.g., a+b | 2*u) and rate constants using the functions k (for forward rate constants) and r (for reverse rate constants), e.g., (a+b | 2*u).k(1.5).r(0.5).

See functions crn_to_odes to convert reactions to ODEs (ordinary differential equations), integrate_crn_odes to get the trajectories of integrating these ODEs over time, and plot_crn to plot the trajectories. The documentation for crn_to_odes explains how reactions are converted into ODEs by each of these functions.

Also supported are inhibitors, which can be added to reactions using the method Reaction.i:

a, b, u, i = gpac.species('A B U I')
rxn = (a+b | 2*u).i(i, 100)

which represents the reaction \(A+B \to 2U\) with inhibitor \(I\) and inhibitor constant 100. Currently the inhibitor is modeled using a first-order Hill function, i.e., its contribution to the reaction rate is to divide by \(1 + i \cdot I\), where \(i\) is the inhibitor constant. So for the reaction defined above, its rate is \([A] \cdot [B] / (1 + 100 \cdot [I])\).

There is also support for simulating discrete chemical reaction networks using the Gillespie algorithm. See the functions plot_gillespie, rebop_crn_counts, and rebop_sample_future_configurations.

KeyConfigCrn module-attribute 

KeyConfigCrn = TypeVar('KeyConfigCrn', Specie, 'sympy.Symbol')

A type variable representing the key type of the dictionary used to represent configurations of CRNs. Typically this is a Specie object, but it can also be a sympy Symbol, for instance when using the dependent_symbols parameter in the functions plot_crn and plot_gillespie, since dependent symbols do not represent species, but rather functions of concentrations of some species.

empty module-attribute 

empty = Expression([])

Used for chemical reactions with empty reactant or product lists, e.g., to implement the exponential decay reaction \(X \to \emptyset\): 

x = species('X')
rxn = x >> empty

Specie dataclass

Represents species in a chemical reaction network. In general these are not created directly, but rather via the species function, which creates a tuple of Specie objects.

name instance-attribute 

name: str

Name of the species. This is used in two ways: when plotting, this name is used in the legend, and when using the dependent_symbols parameter in the functions plot_crn and plot_gillespie, this name is used to identify the species, since dependent symbols need to be specified as sympy.Symbol objects defined as functions of other sympy.Symbol objects. The way to connect that to the species is to make two objects, one Specie and one Symbol, with the same name. See the example notebook for an example of how this is done.

Expression dataclass

Class used for very basic symbolic manipulation of left/right hand side of stoichiometric equations. Not very user friendly; users should just use the species function and manipulate Specie objects with operators >>, |, +, and * to create reactions (see Reaction for examples).

species instance-attribute 

species: list[Specie]

ordered list of species in expression, e.g, A+A+B is [A,A,B]

__getitem__ 

__getitem__(idx: int) -> Specie

Parameters:

Name Type Description Default
idx int

index of species to return

 required

Returns:

Type Description
Specie

Specie at index idx in this Expression

__add__ 

__add__(other: Expression | Specie) -> Expression

Parameters:

Name Type Description Default
other Expression | Specie

Expression or Specie to add to this one

 required

Returns:

Type Description
Expression

Expression representing the union of this Expression and other

__rmul__ 

__rmul__(coeff: int) -> Expression

Parameters:

Name Type Description Default
coeff int

coefficient to multiply this Expression by

 required

Returns:

Type Description
Expression

Expression representing this Expression multiplied by coeff

get_species 

get_species() -> set[Specie]

Returns the set of species in this expression, not their coefficients.

species_counts 

species_counts(key_type: Literal['str', 'Specie'] = 'Specie') -> dict[Specie, int]

Returns a dictionary mapping each species in this expression to its coefficient.

Reaction dataclass

Representation of a stoichiometric reaction using a pair of Expression's, one for the reactants and one for the products.

Reactions are constructed by creating objects of type Specie and using the operators >> (for irreversible reactions) and | (for reversible reactions), as well as the + and * operators to specify the stoichiometric coefficients of the reactants and products, and optionally the methods Reaction.k and Reaction.r to specify forward and reverse rate constants.

For example, the following code creates a reaction that represents the irreversible reaction \(A + B \rightarrow C\) (with implicit rate constant 1.0):

a,b,c = species('A B C')
rxn = a+b >> c

To create reactions

\[ \begin{align*} A+B &\underset{4.1}{\stackrel{0.6}{\rightleftharpoons}} 2C \\ C &\xrightarrow{5.2} D \end{align*} \]

use the following code:

a,b,c,d = gpac.species('A B C D')
rxns = [
    (a+b | 2*c).k(0.6).r(4.1),
    (c >> d).k(5.2),
]

Also supported are inhibitors, which can be added to reactions using the method Reaction.i:

a, b, u, i = species('A B U I')
rxn = (a+b | 2*u).i(i, 100)

which represents the reaction \(A+B \to 2U\) with inhibitor \(I\) and inhibitor constant 100. Currently the inhibitor is modeled using a first-order Hill function, i.e., its contribution to the reaction rate is to divide by \(1 + i \cdot I\), where \(i\) is the inhibitor constant. So for the reaction defined above, its rate is \([A] \cdot [B] / (1 + 100 \cdot [I])\).

reactants instance-attribute 

reactants: Expression = reactants

The left side of species in the reaction.

products instance-attribute 

products: Expression = products

The right side of species in the reaction.

rate_constant class-attribute instance-attribute 

rate_constant: float = k

Rate constant of forward reaction.

rate_constant_reverse class-attribute instance-attribute 

rate_constant_reverse: float = r

Rate constant of reverse reaction (only used if Reaction.reversible is true).

reversible class-attribute instance-attribute 

reversible: bool = reversible

Whether reaction is reversible, i.e. products :math:\to reactants is a reaction also.

inhibitors class-attribute instance-attribute 

inhibitors: list[Specie] = []

Inhibitors of the reaction.

__init__ 

__init__(
    reactants: Specie | Expression,
    products: Specie | Expression,
    k: float = 1,
    r: float = 1,
    reversible: bool = False,
) -> None

In general this constructor should not be used directly; instead, use the operators >>, |, +, and * to construct reactions. (See description of Reaction for examples.)

Parameters:

Name Type Description Default
reactants Specie | Expression

left side of species in the reaction

 required
products Specie | Expression

right side of species in the reaction

 required
k float

Rate constant of forward reaction

 1
r float

Rate constant of reverse reaction (only used if Reaction.reversible is true

 1
reversible bool

Whether reaction is reversible

 False

with_inhibitor 

with_inhibitor(inhibitor: Specie, constant: float = 1.0) -> Reaction

Parameters:

Name Type Description Default
inhibitor Specie

The inhibitor species

 required
constant float

The inhibitor constant

 1.0

i 

i(inhibitor: Specie, constant: float = 1.0) -> Reaction

alias for Reaction.with_inhibitor

get_ode 

get_ode(specie: Specie, reverse: bool = False) -> Expr

Parameters:

Name Type Description Default
specie Specie

A Specie that may or may not appear in this Reaction.

 required
reverse bool

Whether to interpret this reaction in reverse, i.e., treat products as reactants and vice versa. Raises exception if the reaction is not reversible.

 False

Returns:

Type Description
Expr

sympy expression for the ODE term for the given Specie. For example, if the reaction is \(A+B \to 2C\), then the ODE for \(A\) is \(-k \cdot A \cdot B\), the ODE for B is \(-k \cdot A \cdot B\), and the ODE for C is \(2 \cdot k \cdot A \cdot B\).

is_unimolecular 

is_unimolecular() -> bool

Returns:

Type Description
bool

true if there is one reactant

is_bimolecular 

is_bimolecular() -> bool

Returns:

Type Description
bool

true if there are two reactants

symmetric 

symmetric() -> bool

Returns:

Type Description
bool

true if there are two reactants that are the same species

symmetric_products 

symmetric_products() -> bool

Returns:

Type Description
bool

true if there are two products that are the same species

num_reactants 

num_reactants() -> int

Returns:

Type Description
int

number of reactants

num_products 

num_products() -> int

Returns:

Type Description
int

number of products

num_inhibitors 

num_inhibitors() -> int

Returns:

Type Description
int

number of inhibitors

is_conservative 

is_conservative() -> bool

Returns:

Type Description
bool

true if number of reactants equals number of products

reactant_if_unimolecular 

reactant_if_unimolecular() -> Specie

Returns:

Type Description
Specie

unique reactant if there is only one

Raises:

Type Description
ValueError

if there are multiple reactants

product_if_unique 

product_if_unique() -> Specie

Returns:

Type Description
Specie

unique product if there is only one

Raises:

Type Description
ValueError

if there are multiple products

reactants_if_bimolecular 

reactants_if_bimolecular() -> tuple[Specie, Specie]

Returns:

Type Description
tuple[Specie, Specie]

pair of reactants if there are exactly two

Raises:

Type Description
ValueError

if there are not exactly two reactants

reactant_names_if_bimolecular 

reactant_names_if_bimolecular() -> tuple[str, str]

Returns:

Type Description
tuple[str, str]

pair of reactant names if there are exactly two

Raises:

Type Description
ValueError

if there are not exactly two reactants

products_if_exactly_two 

products_if_exactly_two() -> tuple[Specie, Specie]

Returns:

Type Description
tuple[Specie, Specie]

pair of products if there are exactly two

Raises:

Type Description
ValueError

if there are not exactly two products

product_names_if_exactly_two 

product_names_if_exactly_two() -> tuple[str, str]

Returns:

Type Description
tuple[str, str]

pair of product names if there are exactly two

Raises:

Type Description
ValueError

if there are not exactly two products

k 

k(coeff: float) -> Reaction

Same as Reaction.f.

Parameters:

Name Type Description Default
coeff float

The new reaction coefficient

 required

f 

f(coeff: float) -> Reaction

Changes the reaction coefficient to coeff and returns self.

This is useful for including the rate constant during the construction of a reaction. For example

x, y, z = species("X Y Z")
rxns = [
    (x + y >> z).f(2.5),
    (z >> x).f(1.5),
    (z >> y).f(0.5)),
]

Note that if this is a reversible reaction, this specifies the forward rate constant.

Parameters:

Name Type Description Default
coeff float

The new (forward) reaction coefficient

 required

r 

r(coeff: float) -> Reaction

Changes the reverse reaction reaction rate constant to coeff and returns self.

This is useful for including the rate constant during the construction of a reaction. For example, the following defines a reversible reaction \(X + Y \rightleftharpoons Z\) with forward rate constant 2.5 and reverse rate constant 1.5.

x, y, z = species("X Y Z")
rxn = (x + y | z).k(2.5).r(1.5)

Parameters:

Name Type Description Default
coeff float

The new reverse reaction rate constant

 required

get_species 

get_species() -> tuple[Specie, ...]

Returns:

Type Description
tuple[Specie, ...]

a tuple with the species present in the reactants, products, and inhibitors, in that order.

species 

species(sp: str | Iterable[str]) -> tuple[Specie, ...]

Create a tuple of Specie (Single species Expression's).

The return type is tuple[Specie, ...], even if only a single species is found. So the return value always needs to be assigned to a tuple of variables, e.g., x, = species('X').

Examples:

w, x, y, z = species('W X Y Z')
rxn = x + y >> z + w

w, x, y, z = species(['W', 'X', 'Y', 'Z'])
rxn = x + y >> z + w

Note

x, = species('X')
rxn = x >> 2*x

Parameters:

Name Type Description Default
sp str | Iterable[str]

A string or Iterable of strings representing the names of the species being created. If a single string, species names are interpreted as space-separated.

 required

Returns:

Type Description
tuple[Specie, ...]

tuple of Specie objects

crn_to_odes 

crn_to_odes(
    rxns: Iterable[Reaction], species_listed_first: Iterable[Specie] = ()
) -> dict[Symbol, Expr]

Given a set of chemical reactions, return the corresponding ODEs.

Each reaction contributes one term to the ODEs for each species produced or consumed in it. The term from a reaction appearing in the ODE for species X is the product of:

  • the rate constant,
  • the reactant concentrations, and
  • the net stoichiometry of species X in the reaction (i.e., the net amount of X produced by the reaction, negative if consumed).

For example, consider the following two reactions with respective rate constants :math:k_1 and :math:k_2:

\[ \begin{align*} X+X &\xrightarrow{k_1} C \\ C+X &\xrightarrow{k_2} C+Y \end{align*} \]

The net stoichiometry of X in the first reaction is -2, since two copies of X are consumed, and the net stoichiometry of C in that reaction is 1, since one copy of C is produced. The net stoichiometry of C in the second reaction is 0, since it is a catalyst (neither produced nor consumed).

This corresponds to ODEs (following the convention of lowercase letter x for the concentration of species X):

\[ \begin{align*} x' &= -2 k_1 x^2 - k_2 c x \\ c' &= k_1 x^2 \\ y' &= k_2 c x \end{align*} \]

In the package, this can be implemented (for example setting :math:k_1 = 1.5 and :math:k_2 = 0.2) via:

x, y, c = species('X Y C')
rxns = [
    (x+x >> c).k(1.5),
    (c+x >> c+y).k(0.2),
]
odes = crn_to_odes(rxns)
for symbol, ode in odes.items():
    print(f"{symbol}' = {ode}")

which prints

X' = -0.2*C*X - 3.0*X**2
C' = 1.5*X**2
Y' = 0.2*C*X

Parameters:

Name Type Description Default
rxns Iterable[Reaction]

See documentation for Reaction's for details on how to specify reactions.

 required
species_listed_first Iterable[Specie]

If not specified, species are put into the ODEs in the order they are encountered in the reactions. For instance, the reactions X+Y --> A+B and A+X --> L+K would order the species as X, Y, A, B, L, K. But if species_listed_first is [K,Y], then the order would be K, Y, X, A, B, i.e., first listing the species in species_listed_first in the order they are given, and then listing the rest of the species in the order they are encountered in the reactions. This is useful in combination with integrate_odes, since the species trajectories in the returned OdeResult are ordered according to the order of the species in the dict of ODEs given to integrate_odes, so for instance if there are only a two species you care about analyzing, you can list them in species_listed_first and then access result.y[0] and result.y[1] to get their trajectories.

 ()

Returns:

Type Description
dict[Symbol, Expr]

dict mapping each species (represented as a sympy Symbol object, rather than a Specie object) to its corresponding ODE (represented as a sympy Expression). This object can be given as the parameter odes to the functions integrate_odes and plot to integrate/plot the ODEs. (this is essentially what the functions integrate_crn_odes and plot_crn do.

integrate_crn_odes 

integrate_crn_odes(
    rxns: Iterable[Reaction],
    inits: ConfigCrn,
    t_eval: Iterable[float] | None = None,
    *,
    t_span: tuple[float, float] | None = None,
    resets: Mapping[Number, ConfigCrn] | None = None,
    method: str | OdeSolver = "RK45",
    dense_output: bool = False,
    events: Callable | Iterable[Callable] | None = None,
    vectorized: bool = False,
    args: tuple | None = None,
    **options
) -> OdeResult

Integrate the ODEs derived from to the given set of chemical reactions. This calls integrate_odes with the ODEs derived from the given reactions via crn_to_odes.

Parameters:

Name Type Description Default
rxns Iterable[Reaction]

list of Reaction's comprising the chemical reaction network. See documentation for Reaction for details on how to specify reactions.

 required
inits ConfigCrn

dict mapping each species to its initial concentration. Note that unlike the parameter inits in integrate_odes, keys in this dict must be Specie objects, not strings or sympy symbols. Also, symbols explicitly listed in inits also influence the order of trajectories in the OdeResult returned by this function. So while you can leave out a species with initial concentration 0, a reason to specify this explicitly is to make it simpler to access that species' trajectory in result.y, e.g., if inits is {a:0, b:1, c:0, d:2}, this specifies that only b and d have positive initial concentration, an all others are 0. But if result is the OdeResult returned by this function, then result.y[0] corresponds to the trajectory of a, result.y[1] to b, result.y[2] to c, result.y[3] to d. The remaining species are listed in the order in which they are encountered in the reactions.

 required
t_eval Iterable[float] | None

See integrate_odes.

 None
t_span tuple[float, float] | None

See integrate_odes.

 None
resets Mapping[Number, ConfigCrn] | None

If specified, this is a dict mapping times to "configurations" (i.e., dict mapping species to concentrations). The configurations are used to set the concentrations of the species manually during the ODE integration at specific times. Any species not appearing as keys in resets are left at their current values. The OdeResult returned (the one returned by solve_ivp in scipy) will have two additional fields: reset_times and reset_indices, which are lists of the times and indices in sol.t corresponding to the times when the resets were applied. Raises a ValueError if any time lies outside the integration interval, or if resets is empty.

 None
method str | OdeSolver

See solve_ivp

 'RK45'
dense_output bool

See solve_ivp

 False
events Callable | Iterable[Callable] | None

See solve_ivp

 None
vectorized bool

See solve_ivp

 False
args tuple | None

See solve_ivp

 None
options

See integrate_odes.

 {}

Returns:

Type Description
OdeResult

The result of the integration. See integrate_odes for details about this parameter.

plot_crn 

plot_crn(
    rxns: Iterable[Reaction],
    inits: ConfigCrn,
    t_eval: Iterable[float] | None = None,
    *,
    t_span: tuple[float, float] | None = None,
    resets: Mapping[Number, ConfigCrn] | None = None,
    dependent_symbols: dict[Symbol, ValOde] | None = None,
    figsize: tuple[float, float] = default_figsize,
    symbols_to_plot: (
        Iterable[Symbol | Specie]
        | Iterable[Sequence[Symbol | Specie]]
        | str
        | Pattern
        | Iterable[Pattern]
        | None
    ) = None,
    show: bool = False,
    legend: dict[Symbol, str] | None = None,
    latex_legend: bool = False,
    omit_legend: bool = False,
    method: str | OdeSolver = "RK45",
    dense_output: bool = False,
    events: Callable | Iterable[Callable] | None = None,
    vectorized: bool = False,
    return_ode_result: bool = False,
    args: tuple | None = None,
    loc: str | tuple[float, float] = "best",
    warn_change_dpi: bool = False,
    **options
) -> OdeResult | None

Plot the ODEs derived from to the given set of chemical reactions. This calls plot with the ODEs derived from the given reactions via crn_to_odes.

See integrate_crn_odes, integrate_odes, and plot for description of parameters shared with those functions. As with plot, the keyword arguments in options are passed to matplotlib.pyplot.plot, as well as to scipy.integrate.solve_ivp, and as with plot, keyword arguments not recognized by solve_ivp (such as those intended for matplotlib.pyplot.plot) cause solve_ivp to print a warning that it does not recognize the argument.

Note that the parameter dependent_symbols should use sympy symbols, not Specie objects. Here is an example of how to use this parameter. Each species that a dependent symbol depends on should be represented by a sympy symbol with the same name as the corresponding Specie object:

Xp,Xm,Yp,Ym = gpac.species('Xp Xm Yp Ym')
x,y,xp,xm,yp,ym = sympy.symbols('x y Xp Xm Yp Ym')

# dual-rail CRN implementation of sine/cosine oscillator
# x' = -y
# y' = x
rxns = [
    Yp >> Yp + Xm,
    Ym >> Ym + Xp,
    Xp >> Xp + Yp,
    Xm >> Xm + Ym,
    Xp+Xm >> gpac.empty,
    Yp+Ym >> gpac.empty,
]
inits = { Xp: 1, Yp: 0 }
from math import pi
t_eval = np.linspace(0, 6*pi, 200)

dependent_symbols = {
    x: xp - xm,
    y: yp - ym,
}

gpac.plot_crn(rxns, inits, t_eval, dependent_symbols=dependent_symbols, symbols_to_plot=[x,y])

Parameters:

Name Type Description Default
rxns Iterable[Reaction]

See integrate_crn_odes.

 required
inits ConfigCrn

See integrate_crn_odes.

 required
t_eval Iterable[float] | None

See integrate_odes.

 None
t_span tuple[float, float] | None

See integrate_odes.

 None
resets Mapping[Number, ConfigCrn] | None

See integrate_crn_odes.

 None
dependent_symbols dict[Symbol, ValOde] | None

dict mapping each sympy Symbol to an expression that defines its value in terms of other symbols. Note that these are not Specie objects as in the parameter rxns, but sympy symbols. Symbols used in the expression values must have the same name as Specie objects in rxns. Symbols used as keys should not have the same name as any Specie objects in rxns. See the example notebook for an example of how to use this parameter.

 None
method str | OdeSolver

See solve_ivp

 'RK45'
dense_output bool

See solve_ivp

 False
events Callable | Iterable[Callable] | None

See solve_ivp

 None
vectorized bool

See solve_ivp

 False
args tuple | None

See solve_ivp

 None
figsize tuple[float, float]

See plot.

 default_figsize
symbols_to_plot Iterable[Symbol | Specie] | Iterable[Sequence[Symbol | Specie]] | str | Pattern | Iterable[Pattern] | None

See plot. However, one can also use Specie objects directly, which are converted to sympy.Symbol objects with the same name before being passed to plot.

 None
legend dict[Symbol, str] | None

See plot.

 None
latex_legend bool

See plot.

 False
omit_legend bool

See plot.

 False
show bool

See plot.

 False
return_ode_result bool

See plot.

 False
loc str | tuple[float, float]

See plot.

 'best'
warn_change_dpi bool

See plot.

 False
options

See integrate_odes.

 {}

Returns:

Type Description
OdeResult | None

None, or the result of the integration, which is the same as the result of integrate_odes if return_ode_result is True. See integrate_odes for details about this parameter.

rebop_crn_counts 

rebop_crn_counts(
    rxns: Iterable[Reaction],
    inits: Mapping[KeyConfigCrn, int],
    tmax: float,
    *,
    nb_steps: int = 0,
    vol: float | None = None,
    resets: Mapping[Number, Mapping[KeyConfigCrn, int]] | None = None,
    dependent_symbols: Mapping[Symbol, ValOde] | None = None,
    seed: int | None = None
) -> Dataset

Run the reactions using the rebop package for discrete simulation using the Gillespie algorithm.

Parameters:

Name Type Description Default
rxns Iterable[Reaction]

list of Reaction's comprising the chemical reaction network. See documentation for Reaction for details on how to specify reactions.

 required
inits Mapping[KeyConfigCrn, int]

dict mapping each species to its initial integer count. Note that unlike the parameter inits in integrate_odes, keys in this dict must be Specie objects, not strings or sympy symbols. instead of fixed, evenly-spaced time points.

 required
tmax float

the maximum time for the simulation.

 required
nb_steps int

Number of evenly-spaced time points at which to record the counts between 0 and tmax. If not specified, all reaction events and their exact times are recorded,

 0
vol float | None

the volume of the system. If not specified, the volume is assumed to be the sum of the initial counts. reactions with k total reactants have their rate divided by vol^(k-1) to account for the volume.

 None
resets Mapping[Number, Mapping[KeyConfigCrn, int]] | None

If specified, this is a dict mapping times to "configurations" (i.e., dict mapping symbols/str to values). The configurations are used to set the values of the symbols manually during the CRN simulation at specific times. Any symbols not appearing as keys in resets are left at their current values. Raises a ValueError if any time lies outside the simulation interval, or if resets is empty.

 None
dependent_symbols Mapping[Symbol, ValOde] | None

See integrate_crn_odes, except that values must be integers, whereas they can be floats in integrate_crn_odes.

 None
seed int | None

If specified, this is the seed used to initialize the random number generator for the Gillespie algorithm. If not specified, a random seed is generated.

 None

Returns:

Type Description
Dataset

Same Result object returned by rebop.Gillespie.run. (an xarray.Dataset object) It can be indexed by species name to get the counts, and by the key "time" to get the times at which the counts were recorded.

plot_gillespie 

plot_gillespie(
    rxns: Iterable[Reaction],
    inits: Mapping[KeyConfigCrn, int],
    tmax: float,
    *,
    nb_steps: int = 0,
    seed: int | None = None,
    resets: Mapping[Number, Mapping[KeyConfigCrn, int]] | None = None,
    dependent_symbols: Mapping[Symbol, ValOde] | None = None,
    figsize: tuple[float, float] = default_figsize,
    latex_legend: bool = False,
    symbols_to_plot: (
        Iterable[Symbol | Specie]
        | Iterable[Sequence[Symbol | Specie]]
        | str
        | Pattern
        | Iterable[Pattern]
        | None
    ) = None,
    legend: dict[Symbol, str] | None = None,
    omit_legend: bool = False,
    show: bool = False,
    return_simulation_result: bool = False,
    loc: str | tuple[float, float] = "best",
    warn_change_dpi: bool = False,
    vol: float | None = None,
    **options: dict[str, object]
) -> Dataset | None

Similar to plot_crn, but uses the rebop package for discrete simulation using the Gillespie algorithm instead of continuous ODEs.

Undocumented arguments have the same meaning as with plot_crn.

Arguments tmax, nb_steps, and seed are passed to rebop_crn_counts. Any custom keyword arguments (specified in **options) to the function matplotlib.pyplot.plot.

Parameters:

Name Type Description Default
rxns Iterable[Reaction]

See rebop_crn_counts.

 required
inits Mapping[KeyConfigCrn, int]

See rebop_crn_counts.

 required
vol float | None

See rebop_crn_counts.

 None
tmax float

See rebop_crn_counts.

 required
nb_steps int

See rebop_crn_counts.

 0
return_simulation_result bool

Whether to return the simulation result; if True, the result of the simulation is returned, but the default behavior is not to do this, so that if the last line of a notebook cell is a call to plot_gillespie, the result is not printed to the output, only the plot is shown.

 False
seed int | None

See rebop_crn_counts.

 None
resets Mapping[Number, Mapping[KeyConfigCrn, int]] | None

See rebop_crn_counts.

 None
dependent_symbols Mapping[Symbol, ValOde] | None

See rebop_crn_counts.

 None
figsize tuple[float, float]

See plot.

 default_figsize
latex_legend bool

See plot.

 False
symbols_to_plot Iterable[Symbol | Specie] | Iterable[Sequence[Symbol | Specie]] | str | Pattern | Iterable[Pattern] | None

See plot_crn.

 None
legend dict[Symbol, str] | None

See plot.

 None
omit_legend bool

See plot.

 False
show bool

See plot.

 False
loc str | tuple[float, float]

See plot.

 'best'
warn_change_dpi bool

See plot.

 False
options dict[str, object]

See integrate_odes.

 {}

Returns:

Type Description
Dataset | None

The result of the simulation, which is the same as the result of rebop_crn_counts.

rebop_sample_future_configurations 

rebop_sample_future_configurations(
    rxns: Iterable[Reaction],
    inits: Mapping[KeyConfigCrn, int],
    tmax: float,
    trials: int,
    *,
    vol: float | None = None,
    resets: Mapping[Number, Mapping[KeyConfigCrn, int]] | None = None,
    dependent_symbols: dict[Symbol, ValOde] | None = None,
    seed: int | None = None
) -> DataFrame

Sample future configurations of the system using the rebop package. Arguments have same meaning as in rebop_crn_counts, except that nb_steps is not used (set to 1), and there is a parameter trials that specifies how many times to sample the future configuration. The returned dictionary maps each species to a 1D numpy array of length trials, representing the sampled configuration counts for that species.

Parameters:

Name Type Description Default
rxns Iterable[Reaction]

See rebop_crn_counts.

 required
inits Mapping[KeyConfigCrn, int]

See rebop_crn_counts.

 required
tmax float

See rebop_crn_counts.

 required
trials int

how many future configurations to sample

 required
vol float | None

See rebop_crn_counts.

 None
resets Mapping[Number, Mapping[KeyConfigCrn, int]] | None

See rebop_crn_counts.

 None
dependent_symbols dict[Symbol, ValOde] | None

See rebop_crn_counts.

 None
seed int | None

This seed is used to generate random seeds for the rebop package, with a different seed used during each trial to give to See rebop_crn_counts. (Otherwise, we'd just be sampling the same run over and over again.) If not specified, a random seed is generated.

 None

Returns:

Type Description
DataFrame

A pandas DataFrame with one column per species name, giving and number of rows equal to trials, representing the sampled configurations.

replace_reversible_rxns 

replace_reversible_rxns(rxns: Iterable[Reaction]) -> list[Reaction]

Parameters:

Name Type Description Default
rxns Iterable[Reaction]

list of Reaction's

 required

Returns:

Type Description
list[Reaction]

list of Reaction's, where every reversible reaction in rxns has been replaced by two irreversible reactions, and all others have been left as they are

format_rxns 

format_rxns(rxns: Iterable[Reaction]) -> str

Formats reactions using more space than for rxn in rxns: print(rxn), using equal width for each reactant and product to help visualize, e.g.,

A1 + B2 + C3  -->  X1
B1 + C2       -->  Y1
A3            -->  Z1 + X1 + Y1

concentration_to_count 

concentration_to_count(concentration: float, volume: float) -> int

Parameters:

Name Type Description Default
concentration float

units of M (molar) = moles / liter

 required
volume float

units of liter

 required

Returns:

Type Description
int

count of molecule with concentration in volume