gpac documentation

gpac.ode

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. See here for a description of GPACs:

It also has support for a very common model governed by polynomial ODEs, the of continuous mass-action chemical reaction networks:

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.

ode.integrate_odes(odes, initial_values, t_eval=None, t_span=None, dependent_symbols=(), method='RK45', dense_output=False, events=None, vectorized=False, args=None, **options)[source]

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

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,
}
initial_values = {
    a: 10,
    b: 1,
    c: 1,
}
t_eval = np.linspace(0, 1, 5)
gpac.integrate_odes(odes, initial_values, 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):

# 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,
}
initial_values = {
    a: 0,
    b: 1,
    t: 0,
}
t_eval = np.linspace(0, 2*pi, 5) # [0, pi/2, pi, 3*pi/2, 2*pi]
solution = gpac.integrate_odes(odes, initial_values, 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:
Returns:

solution to the ODEs, same as object returned by solve_ivp in scipy.integrate https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

Return type:

OdeResult

ode.plot(odes, initial_values, t_eval=None, t_span=None, dependent_symbols=None, figure_size=(10, 3), symbols_to_plot=None, show=False, method='RK45', dense_output=False, events=None, vectorized=False, args=None, loc='best', **options)[source]

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.) See integrate_odes() for description of parameters below that are not documented.

Parameters:

  • figure_size (Tuple[float, float]) – pair (width, height) of the figure

  • symbols_to_plot (Iterable[Symbol | str] | None) – symbols to plot; if not specified, then all symbols are plotted

  • 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.

  • dependent_symbols (Dict[Symbol | str, Expr | str] | 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.

  • loc (str | Tuple[float, float]) – location of the legend; see documentation for matplotlib.pyplot.legend: https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.legend.html

  • 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, initial_values, t_eval=t_eval, rtol=1e-6)

    For solver-specific parameters to solve_ivp, see documentation for solve_ivp in scipy.integrate: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

    Also used for keyword options to plot in matplotlib.pyplot: https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html. However, note that using such arguments here will cause solve_ivp to print a warning that it does not recognize the keyword argument.

  • odes (Dict[Symbol | str, Expr | str | float]) –

  • initial_values (Dict[Symbol | str, float]) –

  • t_eval (Iterable[float] | None) –

  • t_span (Tuple[float, float] | None) –

  • method (str | OdeSolver) –

  • dense_output (bool) –

  • events (Callable | Iterable[Callable] | None) –

  • vectorized (bool) –

  • args (Tuple | None) –

Return type:

None

gpac.crn

Module for expressing chemical reaction networks and deriving their ODEs. Ideas and much code taken from https://github.com/enricozb/python-crn.

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

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

we can write

a, b, u = species('A B U')
rxns = [
    a+b >> 2*u,
    a+u >> 2*a,
    b+u >> 2*b,
]
initial_values = {a: 0.51, b: 0.49}
t_eval = np.linspace(0, 10, 100)
gpac.plot_crn(rxns, initial_values, 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 = 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 + k \cdot I\), where \(k\) is the inhibitor constant. So for the reaction defined above, its rate is \(k \cdot [A] \cdot [B] / (1 + 100 \cdot [I])\).

crn.crn_to_odes(rxns)[source]

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 \(k_1\) and \(k_2\):

\[ \begin{align}\begin{aligned}X+X &\xrightarrow{k_1} C\\C+X &\xrightarrow{k_2} C+Y\end{aligned}\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}\begin{aligned}x' &= -2 k_1 x^2 - k_2 c x\\c' &= k_1 x^2\\y' &= k_2 c x\end{aligned}\end{align} \]

In the package, this can be implemented (for example setting \(k_1 = 1.5\) and \(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:

rxns (Iterable[Reaction]) – list of Reaction’s comprising the chemical reaction network. See documentation for Reaction for details on how to specify reactions.

Returns:

Dictionary 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 ode.integrate_odes() and ode.plot() to integrate/plot the ODEs. (which is essentially all the functions integrate_crn_odes() and plot_crn() do).

Return type:

Dict[Symbol, Expr]

crn.integrate_crn_odes(rxns, initial_values, t_eval=None, t_span=None, method='RK45', dense_output=False, events=None, vectorized=False, args=None, **options)[source]

Integrate the ODEs derived from to the given set of chemical reactions. This calls ode.integrate_odes() with the ODEs derived from the given reactions via crn_to_odes(). See ode.integrate_odes() for description of parameters other than rxns and initial_values.

Parameters:

  • rxns (Iterable[Reaction]) – list of Reaction’s comprising the chemical reaction network. See documentation for Reaction for details on how to specify reactions.

  • initial_values (Dict[Specie, float]) – dict mapping each species to its initial concentration. Note that unlike the parameter initial_values in ode.integrate_odes(), keys in this dict must be Specie objects, not strings or sympy symbols.

  • t_eval (Iterable[float] | None) –

  • t_span (Tuple[float, float] | None) –

  • method (str | OdeSolver) –

  • dense_output (bool) –

  • events (Callable | Iterable[Callable] | None) –

  • vectorized (bool) –

  • args (Tuple | None) –

Return type:

OdeResult

crn.plot_crn(rxns, initial_values, t_eval=None, t_span=None, dependent_symbols=None, figure_size=(10, 3), symbols_to_plot=None, show=False, method='RK45', dense_output=False, events=None, vectorized=False, args=None, loc='best', **options)[source]

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

See crn.integrate_crn_odes(), ode.integrate_odes(), and ode.plot() for description of parameters. As with ode.plot(), the keyword arguments in options are passed to matplotlib.pyplot.plot (https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html), as well as to scipy.integrate.solve_ivp (https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html), and as with ode.plot(), keyword arguments not recognized by scipy.integrate.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:

  • dependent_symbols (Dict[Symbol | str, Expr | str] | None) – dict mapping each 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 expressions must have the same name as Specie objects in rxns.

  • rxns (Iterable[Reaction]) –

  • initial_values (Dict[Specie, float]) –

  • t_eval (Iterable[float] | None) –

  • t_span (Tuple[float, float] | None) –

  • figure_size (Tuple[float, float]) –

  • symbols_to_plot (Iterable[Symbol | str] | None) –

  • show (bool) –

  • method (str | OdeSolver) –

  • dense_output (bool) –

  • events (Callable | Iterable[Callable] | None) –

  • vectorized (bool) –

  • args (Tuple | None) –

  • loc (str | Tuple[float, float]) –

Return type:

None

crn.species(sp)[source]

Create a list of Specie (Single species Expression’s), or a single one.

Parameters:

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.

Return type:

Tuple[Specie, …]

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
crn.replace_reversible_rxns(rxns)[source]
Parameters:

rxns (Iterable[Reaction]) – list of Reaction’s

Returns:

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

Return type:

List[Reaction]

class crn.Specie(name: 'str')[source]
Parameters:

name (str) –

__init__(name)
Parameters:

name (str) –

Return type:

None

class crn.Expression(species)[source]

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).

Parameters:

species (List[Specie]) –

species: List[Specie]

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

get_species()[source]

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

Return type:

Set[Specie]

__init__(species)
Parameters:

species (List[Specie]) –

Return type:

None

crn.empty = Expression(species=[])

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
crn.concentration_to_count(concentration, volume)[source]
Parameters:

  • concentration (float) – units of M (molar) = moles / liter

  • volume (float) – units of liter

Returns:

count of molecule with concentration in volume

Return type:

int

class crn.Reaction(reactants, products, k=1, r=1, reversible=False)[source]

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}\begin{aligned}A+B &\underset{4.1}{\stackrel{0.6}{\rightleftharpoons}} 2C\\C &\xrightarrow{5.2} D\end{aligned}\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 + k \cdot I\), where \(k\) is the inhibitor constant. So for the reaction defined above, its rate is \(k \cdot [A] \cdot [B] / (1 + 100 \cdot [I])\).

Parameters:

  • reactants (Expression) –

  • products (Expression) –

  • k (float) –

  • r (float) –

  • reversible (bool) –

__init__(reactants, products, k=1, r=1, reversible=False)[source]

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

Parameters:

  • reactants (Specie | Expression) – left side of species in the reaction

  • products (Specie | Expression) – right side of species in the reaction

  • k (float) – Rate constant of forward reaction

  • r (float) – Rate constant of reverse reaction (only used if Reaction.reversible is true

  • reversible (bool) – Whether reaction is reversible

Return type:

None

reactants: Expression

The left side of species in the reaction.

products: Expression

The right side of species in the reaction.

rate_constant: float = 1.0

Rate constant of forward reaction.

rate_constant_reverse: float = 1.0

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

reversible: bool = False

Whether reaction is reversible, i.e. products \(\to\) reactants is a reaction also.

inhibitors: List[Specie]

Inhibitors of the reaction.

with_inhibitor(inhibitor, constant=1.0)[source]
Parameters:

  • inhibitor (Specie) – The inhibitor species

  • constant (float) – The inhibitor constant

Return type:

Reaction

i(inhibitor, constant=1.0)[source]

alias for Reaction.with_inhibitor()

Parameters:

  • inhibitor (Specie) –

  • constant (float) –

Return type:

Reaction

get_ode(specie, reverse=False)[source]
Parameters:

  • specie (Specie) – A Specie that may or may not appear in this Reaction.

  • 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.

Returns:

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\).

Return type:

Expr

is_unimolecular()[source]

Returns: true if there is one reactant

Return type:

bool

is_bimolecular()[source]

Returns: true if there are two reactants

Return type:

bool

symmetric()[source]

Returns: true if there are two reactants that are the same species

Return type:

bool

symmetric_products()[source]

Returns: true if there are two products that are the same species

Return type:

bool

num_reactants()[source]

Returns: number of reactants

Return type:

int

num_products()[source]

Returns: number of products

Return type:

int

num_inhibitors()[source]

Returns: number of inhibitors

Return type:

int

is_conservative()[source]

Returns: true if number of reactants equals number of products

Return type:

bool

reactant_if_unimolecular()[source]

Returns: unique reactant if there is only one Raises: ValueError if there are multiple reactants

Return type:

Specie

product_if_unique()[source]

Returns: unique product if there is only one Raises: ValueError if there are multiple products

Return type:

Specie

reactants_if_bimolecular()[source]

Returns: pair of reactants if there are exactly two Raises: ValueError if there are not exactly two reactants

Return type:

Tuple[Specie, Specie]

reactant_names_if_bimolecular()[source]

Returns: pair of reactant names if there are exactly two Raises: ValueError if there are not exactly two reactants

Return type:

Tuple[str, str]

products_if_exactly_two()[source]

Returns: pair of products if there are exactly two Raises: ValueError if there are not exactly two products

Return type:

Tuple[Specie, Specie]

product_names_if_exactly_two()[source]

Returns: pair of product names if there are exactly two Raises: ValueError if there are not exactly two products

Return type:

Tuple[str, str]

k(coeff)[source]

Same as Reaction.f().

Parameters:

coeff (float) – float The new reaction coefficient

Return type:

Reaction

f(coeff)[source]

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:

coeff (float) – float The new (forward) reaction coefficient

Return type:

Reaction

r(coeff)[source]

Changes the reverse reactionn 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:

coeff (float) – float The new reverse reaction rate constant

Return type:

Reaction

get_species()[source]

Return: the set of species present in the reactants, products, and inhibitors, in the order.

Return type:

Tuple[Specie, …]

Indices and tables