2 - More on models and model description
----------------------------------------

This chapter will focus on the "mini-language" used to describe models
in S-timator and also on the model objects created by function
``read_model()``

.. code:: python

    from stimator import read_model
Picking up the last example from the previous chapter: the *open
two-enzyme system*:

.. figure:: 2chem_open.png
   :alt: Example: a two-reaction open chemical system

   Example: a two-reaction open chemical system

Recall that the minimum components of a model declaration are:

-  title
-  reactions
-  parameters
-  init

.. code:: python

    model_description = """
    title An open two-reaction chemical system
    
    inflow: -> A, rate = kin
    r1: A -> B, rate = k1 * A
    r2: B -> C, rate = k2 * B - k3 * C
    outflow: C ->, rate = kout * C
    
    kin = 0.5
    k1 = 0.1
    k2 = 2
    k3 = 1
    kout = 0.2
    
    init: (A = 0, B = 0, C = 0)
    """
    
    m = read_model(model_description)
.. code:: python

    print m

.. parsed-literal::

    An open two-reaction chemical system
    
    Variables: A B C
    
    inflow:
      reagents: []
      products: [('A', 1.0)]
      stoichiometry: ->  A
      rate = kin
    
    r1:
      reagents: [('A', 1.0)]
      products: [('B', 1.0)]
      stoichiometry: A ->  B
      rate = k1 * A
    
    r2:
      reagents: [('B', 1.0)]
      products: [('C', 1.0)]
      stoichiometry: B ->  C
      rate = k2 * B - k3 * C
    
    outflow:
      reagents: [('C', 1.0)]
      products: []
      stoichiometry: C ->
      rate = kout * C
    
    init: (A = 0.0, C = 0.0, B = 0.0)
    
    k3 = 1
    k2 = 2
    k1 = 0.1
    kin = 0.5
    kout = 0.2
    

2.1 - Reactions
~~~~~~~~~~~~~~~

You can iterate over the reactions of a model:

.. code:: python

    for v in m.reactions:
        print (v)

.. parsed-literal::

    inflow:
      reagents: []
      products: [('A', 1.0)]
      stoichiometry: ->  A
      rate = kin
    
    r1:
      reagents: [('A', 1.0)]
      products: [('B', 1.0)]
      stoichiometry: A ->  B
      rate = k1 * A
    
    r2:
      reagents: [('B', 1.0)]
      products: [('C', 1.0)]
      stoichiometry: B ->  C
      rate = k2 * B - k3 * C
    
    outflow:
      reagents: [('C', 1.0)]
      products: []
      stoichiometry: C ->
      rate = kout * C
    
    

Or, just to get the names:

.. code:: python

    for v in m.reactions:
        print (v.name)

.. parsed-literal::

    inflow
    r1
    r2
    outflow
    

A reaction has a lot of attributes:

.. code:: python

    v = m.reactions.r1
    print (v.name)
    print (v.reagents)
    print (v.products)
    print (v.stoichiometry_string)
    print (v.stoichiometry)
    print (v())

.. parsed-literal::

    r1
    [('A', 1.0)]
    [('B', 1.0)]
    A ->  B
    [('A', -1.0), ('B', 1.0)]
    k1 * A
    

``Model.varnames`` is a list of variable names and ``Model.parameters``
can be used to iterate over the parameters of a model.

.. code:: python

    print (m.varnames)
    for p in m.parameters:
        print ('%6s = %g' % (p.name, p))

.. parsed-literal::

    ['A', 'B', 'C']
        k3 = 1
        k2 = 2
        k1 = 0.1
       kin = 0.5
      kout = 0.2
    

2.2 - Transformations
~~~~~~~~~~~~~~~~~~~~~

Transformations are declared starting a line with a ``~``. These are
quantities that vary over time but are not decribed by differential
equations. In this example ``total`` is a transformation.

.. code:: python

    model_description = """
    title An open two-reaction chemical system
    
    inflow: -> A, rate = kin
    r1: A -> B, rate = k1 * A
    r2: B -> C, rate = k2 * B - k3 * C
    outflow: C ->, rate = kout * C
    
    kin = 0.5
    k1 = 0.1
    k2 = 2
    k3 = 1
    kout = 0.2
    
    init: (A = 0, B = 0, C = 0)
    
    ~ total = A + B + C
    """
    
    m = read_model(model_description)
    print m

.. parsed-literal::

    An open two-reaction chemical system
    
    Variables: A B C
    
    inflow:
      reagents: []
      products: [('A', 1.0)]
      stoichiometry: ->  A
      rate = kin
    
    r1:
      reagents: [('A', 1.0)]
      products: [('B', 1.0)]
      stoichiometry: A ->  B
      rate = k1 * A
    
    r2:
      reagents: [('B', 1.0)]
      products: [('C', 1.0)]
      stoichiometry: B ->  C
      rate = k2 * B - k3 * C
    
    outflow:
      reagents: [('C', 1.0)]
      products: []
      stoichiometry: C ->
      rate = kout * C
    
    total:
      rate = A + B + C
    
    init: (A = 0.0, C = 0.0, B = 0.0)
    
    k3 = 1
    k2 = 2
    k1 = 0.1
    kin = 0.5
    kout = 0.2
    

.. code:: python

    m.solve(tf=50.0, outputs=["total", 'A', 'B', 'C']).plot(show=True)


.. image:: models_files/models_19_0.png


2.3 - Local parameters in processes
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Parameters can also "belong", that is, being local, to processes.

In this example, both ``r1`` and ``r2`` have local parameters

Notice how thes paraemters are listed and refered to in ``print m``:

.. code:: python

    model_description = """
    title An open two-reaction chemical system
    
    inflow: -> A, rate = kin
    r1: A -> B, rate = k * A, k = 0.1
    r2: B -> C, rate = kf * B - kr * C, kf = 2, kr = 1
    outflow: C ->, rate = kout * C
    
    kin = 0.5
    kout = 0.2
    
    init: (A = 0, B = 0, C = 0)
    
    ~ total = A + B + C
    """
    
    m = read_model(model_description)
    print m

.. parsed-literal::

    An open two-reaction chemical system
    
    Variables: A B C
    
    inflow:
      reagents: []
      products: [('A', 1.0)]
      stoichiometry: ->  A
      rate = kin
    
    r1:
      reagents: [('A', 1.0)]
      products: [('B', 1.0)]
      stoichiometry: A ->  B
      rate = k * A
      Parameters:
        k = 0.1
    
    r2:
      reagents: [('B', 1.0)]
      products: [('C', 1.0)]
      stoichiometry: B ->  C
      rate = kf * B - kr * C
      Parameters:
        kr = 1
        kf = 2
    
    outflow:
      reagents: [('C', 1.0)]
      products: []
      stoichiometry: C ->
      rate = kout * C
    
    total:
      rate = A + B + C
    
    init: (A = 0.0, C = 0.0, B = 0.0)
    
    kin = 0.5
    kout = 0.2
    

But this model is exactly the same has the previous model. The
parameters were just made local. (``plot()`` produces the same results).

.. code:: python

    m.solve(tf=50.0, outputs=["total", 'A', 'B', 'C']).plot()


.. image:: models_files/models_24_0.png


The iteration of the parameters is now different:

.. code:: python

    for p in m.parameters:
        print p.name, p

.. parsed-literal::

    kin 0.5
    kout 0.2
    r1.k 0.1
    r2.kr 1.0
    r2.kf 2.0
    

2.4 - External variables
~~~~~~~~~~~~~~~~~~~~~~~~

An external variable is a parameter that appears in the stoichiometry of
a reaction. It is treated as a constant.

In this example, ``D`` is an external variable:

.. code:: python

    m = read_model("""
    title An open two-reaction chemical system
    
    inflow: D -> A, rate = kin * D
    r1: A -> B, rate = k * A, k = 0.1
    r2: B -> C, rate = kf * B - kr * C, kf = 2, kr = 1
    outflow: C -> E, rate = kout * C
    
    D = 1
    kin = 0.5
    kout = 0.2
    E = 2
    
    init: (A = 0, B = 0, C = 0)
    
    ~ total = A + B + C
    """)
    print m

.. parsed-literal::

    An open two-reaction chemical system
    
    Variables: A B C
    
    External variables: D E
    
    inflow:
      reagents: [('D', 1.0)]
      products: [('A', 1.0)]
      stoichiometry: D ->  A
      rate = kin * D
    
    r1:
      reagents: [('A', 1.0)]
      products: [('B', 1.0)]
      stoichiometry: A ->  B
      rate = k * A
      Parameters:
        k = 0.1
    
    r2:
      reagents: [('B', 1.0)]
      products: [('C', 1.0)]
      stoichiometry: B ->  C
      rate = kf * B - kr * C
      Parameters:
        kr = 1
        kf = 2
    
    outflow:
      reagents: [('C', 1.0)]
      products: [('E', 1.0)]
      stoichiometry: C ->  E
      rate = kout * C
    
    total:
      rate = A + B + C
    
    init: (A = 0.0, C = 0.0, B = 0.0)
    
    E = 2
    kout = 0.2
    D = 1
    kin = 0.5
    

.. code:: python

    m.solve(tf=50.0, outputs=['A', 'B', 'C', 'D', 'E']).plot()


.. image:: models_files/models_30_0.png


2.5 - Declaration of outputs
~~~~~~~~~~~~~~~~~~~~~~~~~~~~

You can use ``!!`` to specify what should go into the solution of the
model:

.. code:: python

    m = read_model("""
    title An open two-reaction chemical system
    
    inflow: D -> A, rate = kin * D
    r1: A -> B, rate = k * A, k = 0.1
    r2: B -> C, rate = kf * B - kr * C, kf = 2, kr = 1
    outflow: C -> E, rate = kout * C
    
    D = 1
    kin = 0.5
    kout = 0.2
    E = 2
    
    init: (A = 0, B = 0, C = 0)
    
    ~ total = A + B + C
    
    !! C D E
    """)
    m.solve(tf=50.0).plot()


.. image:: models_files/models_33_0.png


Or use the ``outputs`` argument to the ``solve()`` function (in the form
of a list of desired outputs):

.. code:: python

    m.solve(tf=50.0, outputs=['total', 'A']).plot()


.. image:: models_files/models_35_0.png


``->`` can be used to specify the values of all the rates of all the
processes.

.. code:: python

    m.solve(tf=50.0, outputs=['->', 'D']).plot()


.. image:: models_files/models_37_0.png


2.6 - Explicit differential equations
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code:: python

    m = read_model("""
    title mass on a spring, frictionless
    
    # F = m * a = m * v' = - k * x
    # by Hooke's law and Newton's law of motion
    
    v' = -(k * x) / m
    x' = v
    
    m = 0.5
    k = 1
    
    init: x = 1
    """)
    m.solve(tf=10.0).plot()


.. image:: models_files/models_39_0.png


.. code:: python

    m = read_model("""
    title mass on a spring, with friction
    
    # using Hooke's law and friction proportional to speed,
    # F = m * a = m * v' = - k * x - b * v
    
    v' = (-k*x - b*v) / m
    x' = v
    
    m = 0.5
    k = 1
    b = 0.5
    
    init: x = 1
    """)
    m.solve(tf=10.0).plot()


.. image:: models_files/models_40_0.png


2.7 - Forcing functions
~~~~~~~~~~~~~~~~~~~~~~~

.. code:: python

    m = read_model("""
    title An open two-reaction chemical system
    
    inflow: D -> A, rate = kin * D * step(t, 10, 1)
    r1: A -> B, rate = k * A, k = 0.1
    r2: B -> C, rate = kf * B - kr * C, kf = 2, kr = 1
    outflow: C -> E, rate = kout * C
    
    D = 1
    kin = 0.5
    kout = 0.2
    E = 2
    
    init: (A = 0, B = 0, C = 0)
    
    !! inflow A B C E
    """)
    m.solve(tf=50.0).plot()


.. image:: models_files/models_42_0.png