Combinations

class antimony_combinations.Combinations(mutually_exclusive_reactions: List[Tuple[AnyStr]] = [], directory: Optional[str] = None)

Builds combinations of SBML model using antimony

Create every combination of core hypothesis and extension hypotheses and creates SBML models using antimony from the tellurium package.

Combinations is designed to be subclassed. The necessary user input is given by overriding core functions and providing hypothesis extensions.

The following methods must be implemented (see below for an example):

  • core__reactions()
  • core__parameters()
  • core__variables()

However the following methods are optional:

  • core__functions()
  • core__events()
  • core__units()

Each of these methods should return a valid antimony string, since these strings are used to build up a full antimony model.

Extension hypotheses are added by adding methods to your subclass that begin with extension_hypothesis__. Any method that begins with extension_hypothesis__ will be picked up and used to combinatorially build sbml models.

Any extension_hypothesis__ method should return an instance of the HypothesisExtension class, which is merely a container for some needed information.

Note

Notice the double underscore after extension_hypothesis

Extension Hypotheses can operate in either additive or replace mode, depending on how the models should be combined. additive is simpler. An extension hypothesis is additive when your reaction doesn’t override another, or make another reaction superflous. Examples of such instances might be when adding a mass action reaction to a preexisting set of mass action reactions.

replace mode on the other hand should be used when your reaction should be used instead of another reaction.

Examples:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
class MyCombModel(Combinations):

    # no __init__ is necessary as we use the __init__ from parent class

    def core__functions(self):
        return ''' '''

    def core__variables(self):
        return '''
        compartment Cell;
        var A in Cell;
        var pA in Cell;
        var B in Cell;
        var pB in Cell;
        var C in Cell;
        var pC in Cell;

        const S in Cell
        '''

    def core__reactions(self):
        return '''
        R1f: A -> pA; k1f*A*S;
        R2f: B -> pB; k2f*B*A;
        R3f: C -> pC; k3f*C*B;
        '''

    def core__parameters(self):
        return '''
        k1f    = 0.1;
        k2f    = 0.1;
        k3f    = 0.1;

        k2b    = 0.1;
        k3b    = 0.1;
        VmaxB  = 0.1;
        kmB    = 0.1;
        VmaxA  = 0.1;
        kmA    = 0.1;
        k4     = 0.1;

        S = 1;
        A = 10;
        pA = 0;
        B = 10;
        pB = 0;
        C = 10;
        pC = 0;
        Cell = 1;
        '''

    def core__units(self):
        return None  # Not needed for now

    def core__events(self):
        return None  # No events needed

    def extension_hypothesis__additive1(self):
        return HypothesisExtension(
            name='AdditiveReaction1',
            reaction='pB -> B',
            rate_law='k2b * pB',
            mode='additive',
            to_replace=None,  # not needed for additive mode
        )

    def extension_hypothesis__additive2(self):
        return HypothesisExtension(
            name='AdditiveReaction2',
            reaction='pC -> C',
            rate_law='k3b * C',
            mode='additive',
            to_replace=None,  # not needed for additive mode
        )

    def extension_hypothesis__replace_reaction(self):
        return HypothesisExtension(
            name='ReplaceReaction',
            reaction='pB -> B',
            rate_law='VmaxB * pB / (kmB + pB)',
            mode='replace',
            to_replace='R2f',  # name of reaction we want to replace
        )

    def extension_hypothesis__feedback1(self):
        return HypothesisExtension(
            name='Feedback1',
            reaction='pA -> A',
            rate_law='VmaxA * pA / (kmA + pA)',
            mode='additive',
            to_replace=None,  # name of reaction we want to replace
        )

    def extension_hypothesis__feedback2(self):
        return HypothesisExtension(
            name='Feedback2',
            reaction='pA -> A',
            rate_law='k4 * pA',  # mass action variant
            mode='additive',
            to_replace=None,  # name of reaction we want to replace
        )

Now that we have built a Combinations subclass we can use it as follows:

>>> project_root = os.path.dirname(__file__)
>>> c = MyCombModel(mutually_exclusive_reactions=[
>>>         ('Feedback1', 'Feedback2')
>>>     ], directory=project_root       # optionally specify project root
>>> )

MyCombModel behaves like an iterator, though it doesn’t store all model topologies on the outset but builds models of the fly as the topology attribute is incremented. Topology always starts on model 0, the core model that doesn’t have additional hypothesis extensions.

>>> print(c)
MyCombModel(topology=0)

The complete set of model topologies is enumerated by the topology attribute. The __len__ method is set to the size of this set, accounting for mutually exclusive topologies, which is a mechanism for reducing the topology space.

>>> print(len(c))
24

You can pick out any of these topologies using the selection operator

>>> print(c[4])
MyCombModel(topology=4)

To see which topologies correspond to which hypothesis extensions we can use antimony_combinations.list_topologies(), which returns a pandas.DataFrame.

>>> c.list_topolgies()
                                                      Topology
ModelID
0                                                     Null
1                                                additive1
2                                                additive2
3                                                feedback1
4                                                feedback2
5                                         replace_reaction
6                                     additive1__additive2
7                                     additive1__feedback1
8                                     additive1__feedback2
9                              additive1__replace_reaction
10                                    additive2__feedback1
11                                    additive2__feedback2
12                             additive2__replace_reaction
13                             feedback1__replace_reaction
14                             feedback2__replace_reaction
15                         additive1__additive2__feedback1
16                         additive1__additive2__feedback2
17                  additive1__additive2__replace_reaction
18                  additive1__feedback1__replace_reaction
19                  additive1__feedback2__replace_reaction
20                  additive2__feedback1__replace_reaction
21                  additive2__feedback2__replace_reaction
22       additive1__additive2__feedback1__replace_reaction
23       additive1__additive2__feedback2__replace_reaction

You can extract all topologies into a list using the antimony_combinations.Combinations.to_list() method.

>>> print(c.to_list()[:4])
[MyCombModel(topology=0),
 MyCombModel(topology=1),
 MyCombModel(topology=2),
 MyCombModel(topology=3)]

You can iterate over the set of topologies

>>> for i in c[:3]:
>>> ... print(i)
MyCombModel(topology=0)
MyCombModel(topology=1)
MyCombModel(topology=2)

Or use the items method, which is similar to dict.items().

>>> for i, model in c.items()[:3]:
>>> ... print(i, model)
0 MyCombModel(topology=0)
1 MyCombModel(topology=1)
2 MyCombModel(topology=2)

Selecting a single model, we can create an antimony string

>>> first_model = c[0]
>>> print(first_model.to_antimony())
model MyCombModelTopology0
    compartment Cell;
    var A in Cell;
    var pA in Cell;
    var B in Cell;
    var pB in Cell;
    var C in Cell;
    var pC in Cell;
    const S in Cell
    R1f: A -> pA; k1f*A*S;
    R2f: B -> pB; k2f*B*A;
    R3f: C -> pC; k3f*C*B;
    k1f = 0.1;
    k2f = 0.1;
    k3f = 0.1;
    S = 1;
    A = 10;
    pA = 0;
    B = 10;
    pB = 0;
    C = 10;
    pC = 0;
    Cell = 1;
end

or a tellurium model

>>> rr = first_model.to_tellurium()
>>> print(rr)
<roadrunner.RoadRunner() {
'this' : 0x555a52c8cb90
'modelLoaded' : true
'modelName' :
'libSBMLVersion' : LibSBML Version: 5.17.2
'jacobianStepSize' : 1e-05
'conservedMoietyAnalysis' : false
'simulateOptions' :
< roadrunner.SimulateOptions()
{
'this' : 0x555a5309cd00,
'reset' : 0,
'structuredResult' : 0,
'copyResult' : 1,
'steps' : 50,
'start' : 0,
'duration' : 5
}>,
'integrator' :
< roadrunner.Integrator() >
  name: cvode
  settings:
      relative_tolerance: 0.000001
      absolute_tolerance: 0.000000000001
                   stiff: true
       maximum_bdf_order: 5
     maximum_adams_order: 12
       maximum_num_steps: 20000
       maximum_time_step: 0
       minimum_time_step: 0
       initial_time_step: 0
          multiple_steps: false
      variable_step_size: false

}>

>>> print(rr.simulate(0, 10, 11))
    time,     [A],     [pA],       [B],    [pB],     [C],    [pC]
 [[    0,      10,        0,        10,       0,      10,       0],
  [    1, 9.04837, 0.951626,   3.86113, 6.13887, 5.27257, 4.72743],
  [    2, 8.18731,  1.81269,   1.63214, 8.36786, 4.07751, 5.92249],
  [    3, 7.40818,  2.59182,  0.748842, 9.25116, 3.64313, 6.35687],
  [    4,  6.7032,   3.2968,  0.370018, 9.62998, 3.45361, 6.54639],
  [    5, 6.06531,  3.93469,  0.195519, 9.80448,  3.3609,  6.6391],
  [    6, 5.48812,  4.51188,  0.109779, 9.89022, 3.31158, 6.68842],
  [    7, 4.96585,  5.03415, 0.0651185, 9.93488,  3.2835,  6.7165],
  [    8, 4.49329,  5.50671, 0.0405951,  9.9594, 3.26657, 6.73343],
  [    9,  4.0657,   5.9343, 0.0264712, 9.97353, 3.25584, 6.74416],
  [   10, 3.67879,  6.32121, 0.0179781, 9.98202, 3.24872, 6.75128]]

Or an interface to copasi, via pycotools3

>>> c.to_copasi()
Model(name=NoName, time_unit=s, volume_unit=l, quantity_unit=mol)

Which could be used to configure parameter estimations. Currently, support for parameter estimation configuration has in COPASI not been included but this is planned for the near future.