{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# MuMoT User Manual \n", "\n", "## Multiscale Modelling Tool \n", "\n", "This is the user manual for [MuMoT](https://mumot.readthedocs.io/) (Marshall *et al.*, [2019](#references)), a software tool developed at the University of Sheffield as part of the [DiODe](http://diode.group.shef.ac.uk) project." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Working with MuMoT\n", "\n", "MuMoT runs inside [Jupyter notebooks](http://jupyter.org) - since you are reading this User Manual you have presumably already installed Jupyter, or are using an installation provided to you. To toggle the [table of contents](https://github.com/ipython-contrib/jupyter_contrib_nbextensions/tree/master/src/jupyter_contrib_nbextensions/nbextensions/toc2) click the 'Table of Contents' button in the Jupyter toolbar.\n", "\n", "Next, you need to install MuMoT itself; MuMoT is a Python package - see the [project documentation](https://mumot.readthedocs.io/en/latest/getting_started.html) for detailed installation instructions.\n", "\n", "To run commands in 'Code' cells it's usually easiest just to hit `Shift + Enter`; here is some basic information on [working with Jupyter notebooks](http://jupyter-notebook.readthedocs.io/en/stable/examples/Notebook/Notebook%20Basics.html). You can also use [Markdown](https://en.wikipedia.org/wiki/Markdown) cells (like this one) to write descriptive text, using simple [formatting](https://www.markdownguide.org/cheat-sheet) and even LaTeX [maths](https://en.wikibooks.org/wiki/LaTeX/Mathematics). Double-click in a Markdown cell to see how it is coded; hit `Shift + Enter` to return to the reader's-eye view of it. There are other kinds of cell too; to change the current cell's type click on the drop-down widget in the Jupyter toolbar.\n", "\n", "With Jupyter and MuMoT both installed, you simply need to import the MuMoT package into your notebook using `import mumot`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import mumot\n", "\n", "mumot.about()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Defining your first model\n", "\n", "If everything went well with importing MuMoT you should have seen the message\n", "\n", "``Created `%%model` as an alias for `%%latex``\n", "\n", "followed by information on the current version, and documentation, from the `about()` command.\n", "\n", "MuMoT models are defined using the `%%model` keyword within a cell, and have a very simple syntax; since MuMoT works with *chemical kinetic* or *reaction kinetic* models, you simply have to describe the reactions that take place, and the rates with which they occur. For example\n", "\n", "`A + A -> A + U: s`\n", "\n", "is a very simple reaction, in which two particles of type $A$ interact, and one changes to type $U$, at rate $s$.\n", "\n", "These kind of reaction rules are quite general; they describe a lot of chemical reactions, as well as collective behaviour in which individuals interact to change each other's state, and even demographic models in which individuals are born and die.\n", "\n", "Our first model will be based on signalling behaviour observed in honeybee swarms, as described and modelled in Seeley *et al.* ([2012](#references)) and Pais *et al.* ([2013](#references)).\n", "\n", "\n", "The `$` at the start and end just allow us to use [LaTeX](https://en.wikipedia.org/wiki/LaTeX) codes for our state and rate labels, allowing us to use letters from different alphabets, such as Greek, and other nice formatting - if you don't know LaTeX don't worry about this - leave in the `$` signs, they won't hurt.\n", "\n", "OK, here we go..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%model\n", "$\n", "U -> A : g_A\n", "U -> B : g_B\n", "A -> U : a_A\n", "B -> U : a_B\n", "A + U -> A + A : r_A\n", "B + U -> B + B : r_B\n", "A + B -> A + U : s\n", "A + B -> B + U : s\n", "$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we have our model, and Jupyter should have reprinted the rules underneath the cell you input them in.\n", "\n", "You can see that actually the `_` code enables you to use subscripts, which can be handy for distinguishing related rates - in this model there are two distinct nest sites for the honeybee swarm to choose between, *A* and *B*, and rates associated with each one.\n", "\n", "We have our model, but currently it is not useful; it's only a definition, not a fully-fledged model that we can start analysing with MuMoT. To make one of those we have to use the `parseModel` command; we simply pass in a reference to the cell where we defined our model - this should be in input cell 2..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model1 = mumot.parseModel(In[2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If there is no output from `parseModel` that is good news - it means the model was valid. Now we have an object, which we called `model1`, which is a version of the model that we can start doing useful work with. For example, we can see its 'reactants' (or states, or types of individual; remember the terminology comes from chemical reactions)..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model1.showReactants()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "...and the rates the reactions happen at..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model1.showRates()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that here MuMoT also helpfully tells us which reactions in the original model definition each rate is associated with.\n", "\n", "We can also look at the reaction equations in a slightly nicer format..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model1.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can even see a simple figure representing our model..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model1.visualise()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This figure attempts to represent the nature of the interactions between reactants graphically, from simple transitions (arrows), inhibitions (filled circles) and induced switches (arrows with filled circles at the origin); to understand the different interaction patterns try relating the arrows with the reactions we defined above, by reference to the rates they are labelled with." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exploring your new model\n", "\n", "If you want to know what commands you can send to an object, Jupyter lets you type out the object name, then a `.`, then hit `Tab` to see a list of accepted commands. You can also consult the [Command reference](#comreference) for a full list of model commands, and details on their use.\n", "\n", "One commmand that looks intruiguing is `showODEs()`" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model1.showODEs()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is where things start getting interesting, because equations are what mathematicians analyse. `showODEs` shows the model's Ordinary Differential Equations, which describe how the populations of reactants in the model change over time. MuMoT has derived these equations for you automatically, from your description of the reactions. If you are familiar with ODEs then you should be able to read these quite easily. However even if you are not, you can still work with them. By the way, if you right click on the equations you can get the source code to generate them in a number of different formats, should you want to insert them into a paper or report.\n", "\n", "So, what analysis can we do in MuMoT now we've found our ODEs?\n", "\n", "Well, we might want to look at how the state of the system (the proportions of the different reactants) change on average over time...\n", "\n", "...but before we do that, we let MuMoT know that this is actually a simpler set of equations than it appears to be. $A$, $B$, and $U$ represent states that honebees can be in during the swarming process, either committed to site *A*, committed to site *B*, or *U*ncommitted. Because there is a constant number of individuals in the swarm during a decision, one of these equations is redundant; the change in the uncommitted bee frequency can be worked out from the change in frequencies of $A$-committed and $B$-committed bees, for example.\n", "\n", "We let MuMoT know this by defining one of the reactants in terms of the other, using `substitute()`" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model2 = model1.substitute('U = N - A - B')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we have a simplified model, `model2`, which only has two equations, because we have told MuMoT that the total number of uncommitted bees $U$ is the total number of bees in the swarm $N$, minus the numbers committed to the two different nests ($A$ and $B$)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model2.showODEs()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Having a closer look at these equations, we see that $U$ has vanished from the equations. \n", "\n", "If you want to find out why the ODEs look the way they do, how they were derived, and where the funny symbols come from, keep reading. If not, skip forward a couple of cells to where we generate out first interactive graphical *view* on the model.\n", "\n", "Let's inspect the ODEs of `model2` (and those of `model1` shown a bit further up) even closer. The symbols denoting *reactants* in the ODEs are *concentrations*, as their derivation is based on the *law of mass action*: the rate with which a chemical reaction occurrs is proportional to (the product of) the concentration of the reactants. This can be made clearer by deriving the ODEs via a principled method based on the system size expansion (also called *van Kampen* expansion ([van Kampen, 2007](#references)) of the *Master equation*, as explained in the [relevant section](#masterEqFokkerPlanck) below. In MuMoT we can use `showODEs()` with the keyword `method='vanKampen'` - here is an example: " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model2.showODEs(method='vanKampen')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, we use the upper-case Greek letter $\\Phi$ to denote concentrations (proportions) and the reactants are added as indices. This is done throughout MuMoT and you will encounter this notation again when you work your way through this notebook. Hence, upper-case $\\Phi$ cannot be used as a user-defined symbol in MuMoT to avoid confusion. \n", "\n", "If the keyword `method='vanKampen'` is not provided in `showODEs()`, a heuristic default method is applied accessible with the keyword `method='massAction'`. As `method='massAction'` is used as default automatically, it does not need to be given as a *keyword* in the `showODEs` command. Using the default method is faster than deriving the ODEs via *van Kampen* expansion. To discriminate between both methods we use the two different notation styles. However, both methods yield equations which describe how the system evolves over time on average.\n", "\n", "This shall be enough about mathematical equations for now; let's turn to more exciting graphical output and numerically *integrate* the ODEs using the MuMoT-method `integrate()`" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "int1 = model2.integrate()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Above you should see an interactive plot of how the reactant populations change over time.\n", "There are various ways to manipulate this graph. You will notice that there are some sliders above the figure; MuMoT gives you a slider for every rate in the model, and by changing their values, you can see what happens to the behaviour of the system. There is also an `Advanced options` tab which unfolds if you click on it. You can choose the initial conditions (note those are given as proportions labelled $\\Phi_{}$), the simulation time and whether you want to plot *absolute numbers* (labelled with the reactants' symbols as given in the original model definition) or *proportions* (labelled $\\Phi_{}$). In this example, the initial conditions you can choose on the sliders sum to 1 automatically; recall that this is because the system is *conserved* (i.e. the total number of bees is a constant). MuMoT recognises this, as one of the original reactants (representing uncommitted bees $U$) was substituted using `substitute()`. There are many more plotting options available which are all explained in the section [Keyword arguments](#sec:kwargs). \n", "\n", "Jupyter also lets you interact nicely with the figure; you can zoom in, back out, pan around, and save, all by clicking the icons that appear underneath it. If you want to save a figure for a paper, simply click the disk icon. And if while exploring the model you find a nice figure, you can click the *off icon* in its top right to snapshot it, and generate a new interactive figure to play with (if you start running out of space for your plots, click in the left margin to toggle scrolling). You can then simply right-click your snapshots to save them from your browser - just don't forget to note the parameters that generated them - you can find that out by looking at the logs... later in the manual we will see how to take notes properly to reproduce analyses (see [Bookmarking](#bookmarking)).\n", "\n", "Every *controller* in MuMoT keeps a log of everything it's been asked to do; we called our controller `int1` when we created it, so it is easy to check what it has been up to... (WARNING: if you have been playing with the sliders a lot, there could be a lot in the logs! To only see the last few entries we set the `tail` argument to `True`; if `tail` is not set it will be assumed to be `False`)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "int1.showLogs(tail = True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another possibility to look at how the model behaves is by asking for a *vector plot* of the two reactants $A$ and $B$ ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "vector1 = model2.vector('A', 'B', showNoise = True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The figure you see is a vector plot, which shows on average how the state of the system (proportions of bees in states $A$ and $B$ denoted $\\Phi_{A}$ and $\\Phi_{B}$ following our convention introduced above) will change over time, from any point in the *state space*. The arrows give the direction the system will move in, and their lengths show how fast. Because the total number of bees is constant, so $A+B$ can never exceed 1, the top right hand triangle of the figure is greyed out; these are impossible states. The number of uncommitted bees $U$, which we hid away earlier, is not shown, but it is simply $1-A-B$. But what we really care about in a decision-making process is how the number of 'voters' for the competing options changes over time.\n", "\n", "Again, model parameters may be altered via the sliders above the figure and useful information about the computation can be obtained using:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "vector1.showLogs(tail = True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For papers and presentations we can produce somewhat prettier plots similar to figure 1 above, using `stream()`. But before we do that, one unfortunate thing about our `vector1` controller is that there are just too many sliders. We might choose to describe some of the rates (*a*bandonment, *r*ecruitment (via the waggle dance), and individual discovery and commitment (*g*)) in terms of the *v*alues of the nest sites... we can do this by using `substitute()` again:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model3 = model2.substitute('a_A = 1/v_A, a_B = 1/v_B, g_A = v_A, g_B = v_B, r_A = v_A, r_B = v_B')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our new model is a little simpler, with just three rates, $v_A$, $v_B$ and $s$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model3.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Still, in decision-making what we really care about is how good options are on average, and their difference in value. Pais *et al.* ([2013](#references)) thus defined option values in terms of deviation from the average value, so\n", "$$\\displaystyle v_A=\\mu+\\frac{\\Delta}{2}$$\n", "and\n", "$$\\displaystyle v_B=\\mu-\\frac{\\Delta}{2}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model4 = model3.substitute('v_A = \\mu + \\Delta/2, v_B = \\mu - \\Delta/2')\n", "model4.showODEs()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This model is much easier to interact with..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "stream1 = model4.stream('A','B', showFixedPoints = True, initWidgets={'mu':[3, 1, 5, 0.5]})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A `stream` plot is just like a `vector` plot, except now lines show the average change of the system over time in more detail, and their shading represents the speed of change, from slow (light grey) to fast (black).\n", "\n", "Note also that here we added a *keyword argument* `showFixedPoints = True` - this can also be applied to `vector` plots, and plots the stationary points of the equation. If the stationary point is *stable* (arrows converge on that point) then it is represented with a filled circle, and if it is *unstable* (arrows leave that point) then a hollow circle is shown. Displaying *fixed points* can slow down plots, so by default they are not shown, and you need to explicitly request they are included. To know if a calculation is still running, the *bookmark* button changes colour. It is *pink* during a computation and *grey* otherwise (maybe you have noticed this already in the plots above). \n", "\n", "Analysis of fixed points is an important technique for understanding *dynamical systems* such as our model. But while exploring a model interactively by changing its parameters and seeing how it behaves can help develop intuition, really we need a more systematic approach. Fortunately we can do this via *bifurcation plots*, using the `bifurcation()` command...\n", "\n", "`bifurcation()` takes two primary arguments, the first is the parameter we want to vary systematically, the second is the state we want to observe. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bifurcation1 = model4.bifurcation('s','A-B', \n", " initWidgets={'mu':[3, 1, 5, 0.5], 'Delta':[0, 0, 2, 0.1], \n", " 'initBifParam':[4.5, 4, 6, 0.1]},\n", " choose_xrange=[0, 5])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our new bifurcation plot lets us see what happens to fixed points in the system as we vary a control parameter, in this case the cross-inhibitory stop-signalling rate between bees, $s$. We have also chosen not to observe a single reactant as our state variable, but rather the difference between two of them; this makes sense because in decision-making we care about the difference in votes for options.\n", "\n", "Note that in the `bifurcation()` command above we specified initial parameter values (and ranges) for the sliders (via `initWidgets={...}`) and the range shown on the x-axis (with `choose_xrange=[xmin, xmax]`), which yields appropriate results straight away. This was necessary because bifurcation methods are very sensitive to even tiny numerical inaccuracies; nice plots are therefore not guaranteed and typically some parameter configurations work better than others. \n", "\n", "Our plot shows the location of fixed points as the *bifurcation parameter* is varied, and they are either depicted as stable (solid line) or unstable (dashed line). Where the stability of fixed points changes, or fixed points appear or dissappear, these are referred to as *bifurcation points*. The label *BP1* means there is one occurrence of a *branch point*. This type of bifurcation is referred to as a *pitchfork bifurcation*. As a next step, move the $\\Delta$-slider to 0.1 and observe an *unfolding* of the *pitchfork bifurcation* where the branches become disconnected and instead of a *branch point* a *limit point* (*LP1*) is now detected (this computation may take some time). Now it's your turn. Using the interactive `bifurcation` plot above see if you can reproduce the bifurcation diagrams shown in figure 5(i) and 5(ii) of Pais *et al.* ([2013](#references)).\n", "\n", "Figure 5(iii) of Pais *et al.* ([2013](#references)) requires us to specify a different bifurcation parameter..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bifurcation2 = model4.bifurcation('\\\\Delta','A-B', \n", " initWidgets={'mu':[3, 1, 5, 0.5], 's':[5, 4, 6, 0.5]},\n", " choose_xrange=[-1,1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "More information about what `bifurcation()` currently can (and cannot) do is given in the relevant [Advanced options](#advancedbif) section." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The effects of finite system size\n", "\n", "The ODE models we have been looking at based on our model are very useful, but they represent an idealised scenario; they describe the behaviour of a system with infinitely many components or, to put it another way, the *expected* or *mean field* behaviour of the system, excluding noise.\n", "\n", "Noise is very important in collective behavour systems, however, and can have a variety of causes. The first cause we will consider is noisy fluctuations from the infinite population ideal, caused by having a finite, much smaller, population. To begin analysing the noise in our system, dependent on its size, we can deploy two main techniques; stochastic simulation, and statistical physics analysis.\n", "\n", "We will return to stochastic simulation later, concentrating first on statistical physics analysis. As you might imagine, the maths involved in a statistic physics analysis is very advanced, based on deriving and then analysing the *Master equation*... it was already mentioned above in discussing how the model's ODEs were derived; fortunately MuMoT automates this analysis for you.\n", "\n", "Brace yourself, this might look intimidating..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model1.showMasterEquation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You really don't need to understand the equation above, but it's nice to know that MuMoT generated it for you, isn't it? The Master equation describes the probabilities of different transitions between system states. We are particularly interested in it because we can use it to derive the infinite-population behaviour of the system (the ODEs we already saw), but also the noise around that. This is done by applying a technique called the *van Kampen expansion* ([van Kampen, 2007](#references)), to approximate the *Fokker-Planck equation*. In the *Advanced users* section at the end of this manual you'll be able to look at these equations, if you dare, to convince yourself that MuMoT is doing some clever stuff for you\n", "\n", "But what is this clever stuff useful for? Well, for one thing, it lets us see approximately what the effects of noise, caused by finite system size, are in terms of deviations from the fixed points of the system..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "stream2 = model4.stream('A','B', showFixedPoints = True, showNoise = True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This new `stream` plot now lets us see how the system's *expected state* at *equilibrium* gets 'spread out' around the fixed points. Slide the system size slider to see how the extent of the spread depends on the number of components in the system; slide the other sliders to see how the shape of the noise changes, in a way we might not have guessed, as the parameters change.\n", "\n", "The effects of noise in the `stream` plot above are based on the computational approximations of the Master equation. When the system has fewer than one or two reactants the noise effects can be approximated mathematically from the Master equation itself; otherwise the computational approximation is achieved through the efficient and accurate *Stochastic Simulation Algorithm* (SSA) (or *Doob-Gillespie algorithm*) [(Gillespie, 1976)](#references). MuMoT also gives you direct access to the SSA simulator to simply visualise the behaviour of the system over time.\n", "\n", "To run the SSA simulations on your model, you can use the command `SSA()`. This command allows you to perform a single simulation run; you can specify several advanced parameters from the *Advanced options* tab, as described in the relevant [Advanced options](#advoptionsSSA) section.\n", "Try modifying the parameters and see how they change your results." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ssa1 = model4.SSA(initWidgets={'initialState': {'U': [1,0,1,0.05], 'B': [0,0,1,0.05], 'A': [0,0,1,0.05]},\n", " 'systemSize':[50,5,100,1] })" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Spatial noise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "MuMoT also enables more sophisticated stochastic simulations of your models. The `multiagent()` command allows you to run simulations either with static agents interacting over a static communication network, or mobile agents communicating to neighbours within a local range of communication.\n", "\n", "Each agent represents a reactant, and exchanges messages to interact with other agents and change its internal state according to a *probabilistic finite state machine*. The conversion of transition rates into probabilities is done following the methodology proposed in Reina *et al.* ([2005](#references)).\n", "\n", "The multiagent simulations have several advanced options, as described in the relevant [Advanced options](#advoptionsmulti) section. Try modifying the parameters to see how they change your results." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "multiagent1 = model4.multiagent(initWidgets={'initialState': {'U': [1,0,1,0.05], 'B': [0,0,1,0.05], 'A': [0,0,1,0.05]},\n", " 'systemSize':[20,5,100,1], 'visualisationType':'graph' })" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Going to a higher dimension\n", "Our first model ended up being 2-dimensional, after we realised that there were a constant number of components in the system. However MuMoT can work with higher-dimensional models (although not all features are supported for them; see [Supported features](#features))\n", "\n", "For example, we can consider the generalisation of `model1` above to 3 or more alternatives, as presented by Reina *et al.* ([2017](#references))...\n", "\n", "(note here that we will use a relative reference to the last active cell, to avoid problems if we add or remove cells to the notebook; due to the peculiarities of Python the last cell has reference `-2`)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%model\n", "$\n", "U -> A : k * v_A\n", "U -> B : k * v_B\n", "U -> C : k * v_B\n", "A -> U : k * (1 / v_A)\n", "B -> U : k * (1 / v_B)\n", "C -> U : k * (1 / v_C)\n", "A + U -> A + A : h * v_A\n", "B + U -> B + B : h * v_B\n", "C + U -> C + C : h * v_C\n", "A + B -> A + U : h * v_A\n", "A + B -> B + U : h * v_B\n", "A + C -> A + U : h * v_A\n", "A + C -> C + U : h * v_C\n", "B + C -> B + U : h * v_B\n", "B + C -> C + U : h * v_C\n", "$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model5 = mumot.parseModel(In[-2]).substitute('U = N - A - B - C')\n", "vector2 = model5.vector('A', 'B', 'C')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Non-constant system size\n", "Not all models are of the kind studied above, in which there are a constant number of components. In *open systems* components can be created and destroyed, so the total size of the system can fluctuate over time, and in principle has no limits.\n", "\n", "A classic model from both ecology and chemistry is the *Lotka-Volterra* model; ecologists often think of this as representing two interacting populations, one of predators, and one of prey. Individuals are born and die in both populations. MuMoT enables this both via the `\\emptyset` notation $\\emptyset$, in which components spontaneously arise from nothing, or are destroyed, and via steady-state reactant populations, denoted with `()`, which do not change over time but simply persist at a constant size.\n", "\n", "Our Lotka-Volterra example is taken from Murray ([2002](#references), p124). Individuals are born into the prey population $X$ at a rate that depends on that population size, then predatory interactions decrease the prey population size and increase the predator population size. Finally, predators experience background mortality." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%model\n", "$\n", "(A) + X -> X + X : \\alpha\n", "X + Y -> Y + Y : \\beta\n", "Y -> \\emptyset : \\gamma\n", "$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model6 = mumot.parseModel(In[-2])\n", "model6.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model6.showODEs()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that there is no equation for $A$; since it is specified to be at a constant concentration, it does not change over time. This means that when we try and interact with the model, we will need to specify the constant value that $A$ should take.\n", "\n", "Also, since the system size is now no longer constant, as discussed above, it is useful to be able to change the limits on the plot, to ensure that all the interesting behaviour of the system is visible... the `Plot limits` widget enables this." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "stream3 = model6.stream('X','Y', showFixedPoints = False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "### Bookmarking\n", "\n", "If you find an interesting parameter set and want to preserve it, you can store those parameter values in a form that lets you reproduce your analysis. Simply click on the bookmark button underneath the sliders; you should see a message that a bookmark has been pasted to the logs, and instructions how to view it. Let's follow these instructions from our latest plot... click the bookmark button on the last plot, then run the following command:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "stream3.showLogs(tail = True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we need simply copy the last line of the log file, paste into a new cell, and replace `` with the model we are working with, in this case `model5`; this is easily done by clicking in the cell then using *Find and Replace* from the *Edit* menu. We will also give our new bookmarked view a unique name, so we can refer to it if we need to in the future... for example" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bookmark1 = model6.stream('X', 'Y', params = [('\\\\beta', 0.8), ('\\\\alpha', 0.5), ('A', 0.5), ('\\\\gamma', 0.5), ('plotLimits', 3.5), ('systemSize', 1)], showFixedPoints = False, bookmark = False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Command reference\n", "\n", "\n", "#### Getting help\n", "You can get help on commands by using the `help()` command. To get basic information on MuMoT, for example, type" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "help(mumot)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To get help on a specific command, such as the `stream` command, available for a model, type" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "help(model1.stream)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Supported features\n", "\n", "##### Graphical output\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Command1d-system2d-system3d-systemempty-setconstant reactants
stream✕✓✕✓✓
vector✕✓✓✓✓
SSA✓✓✓✓✓
multiagent
(static network)
(moving particles)
✓
✓
✓
✓
✓
✓
✕
✓
✕
✓
bifurcation✓✓✕✓✓
integrate✓✓✓✓✓
noiseCorrelations✓✓✓✓✓
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Display of mathematical expressions\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Command1d-system2d-system3d-system4d-systemempty-setconstant reactants
showODEs✓✓✓✓✓✓
showMasterEquation✓✓✓✓✓✓
showVanKampenExpansion✓✓✓✓✓✓
showFokkerPlanckEquation✓✓✓✓✓✓
showNoiseEquations✓✓✓✓✓✓
showNoiseSolutions✓✓✓✓✓✓
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Model properties\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Command1d-system2d-system3d-system4d-systemempty-setconstant reactants
show✓✓✓✓✓✓
showRates✓✓✓✓✓✓
showReactants✓✓✓✓✓✓
showConstantReactants✓✓✓✓✓✓
showStoichiometry✓✓✓✓✓✓
visualise✓✓✓✓✓✓
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Result download\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Command1d-system2d-system3d-system4d-systemempty-setconstant reactants
downloadResults✓✓✓✓✓✓
\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Commandbifurcationintegratemultiagentnoise correlationsSSAstreamvector
downloadResults✕✕✓✕✓✕✕
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "##### Access to symbolic mathematical expressions\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Command1d-system2d-system3d-system4d-systemempty-setconstant reactants
getODEs✓✓✓✓✓✓
getStoichiometry✓✓✓✓✓✓
getMasterEquation✓✓✓✓✓✓
getVanKampenExpansion✓✓✓✓✓✓
getFokkerPlanckEquation✓✓✓✓✓✓
getNoiseEquations✓✓✓✓✓✓
getNoiseSolutions✓✓✓✓✓✓
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Advanced options\n", "\n", "\n", "All advanced options can be set via [keywords](#sec:kwargs)\n", "\n", "All controller-generating commands have the optional arguments `params`, which contains a list of fixed parameter values (see [Partial controllers](#partialcont)), and `initWidgets`, which contains a list of widget configurations (see [Initialising widgets](#initwidgets))\n", "\n", "In addition, all plotting functions accept a number of formatting keywords, those are explained below in the [Keyword arguments](#sec:kwargs) section separately for each plotting function." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "##### Advanced options for `stream()` and `vector()`\n", "* If keyword `showNoise = True` given:\n", " * `maxTime`: simulation time: the length of each simulation run\n", " * `randomSeed`: a number to control (and repeat) random stochastic effects. Any run that uses the \n", " * `runs`: number of simulation runs to be executed\n", " * `aggregateResults`: flag to aggregate the results from several runs\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "##### Advanced options for `bifurcation()`\n", "\n", "* `initBifParam`: starting value of the bifurcation parameter. MuMoT tries to determine all stationary states and starts numerical continuation (using PyDSTool) from all stable fixed points found.\n", "* `initialState`: this will only have an effect if no stationary solutions could be detected. Then these initial conditions will be passed on to PyDSTool for the numerical continuation.\n", "\n", "Note: Currently only 1D and 2D ODE systems are supported in `bifurcation()` and only Branch Points (BP) and Limit Points (LP) can be detected in 1-parameter bifurcations (or more appropriately, codimension one bifurcations)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Advanced options for `SSA` and `multiagent`\n", "\n", "\n", "* `initialState`: proportion of each reactant at timestep zero: remember that the sum of all (non-constant) reactants must be 1.0\n", "* `maxTime`: simulation time: the length of each simulation run\n", "* `randomSeed`: a number to control (and repeat) random stochastic effects. Any run that uses the same random seed will have tha same identical random fluctuations\n", "* `plotProportions`: flag to plot proportions (i.e. reactants are in closed interval [0,1]) or full populations (i.e. reactants are in closed interval [0,systemSize])\n", "* `realtimePlot`: flag to plot results in realtime. If `True`, the plot is updated each timestep of the simulation; if `False`, the plot is updated once at the end of the simulation\n", "* `runs`: number of simulation runs to be executed\n", "* `aggregateResults`: flag to aggregate the results from several runs\n", "* `visualisationType`: type of visualisation:\n", " * `evo`: Temporal evolution with y-axis showing the population size of each reactant varying as a function of the time on the x-axis. When aggregated, the results are visualised as boxplots. \n", " * `final`: Final distribution of the population among the reactants. In this visualisation, the use can specify the reactants on x and y axes. When aggregated, the results are visualised as an ellipse estimated as the covariance matrix\n", " * `barplot`: Barplot showing the final distribution. When aggregated, the barplots show the average and the error bars show the standard deviation; otherwise the results of the last simulation run are displayed\n", "* `final_x`: reactant to plot as X-axis in final distribution visualisation\n", "* `final_y`: reactant to plot as Y-axis in final distribution visualisation\n", "\n", "(Note that modifying the widgets: type of visualisation, flag to plot proportions, and flag to aggregate does not trigger the computation of new data. All other widgets do.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Advanced options for `multiagent` only\n", "\n", "\n", "`multiagent` views have all the options described above for the `SSA` simulations and, in addition:\n", "* `netType`: type of network. Available types are:\n", " * `full`: Full graph; a complete network in which every agent communicates with every agent\n", " * `erdos-renyi`: Random network with Erdos-Renyi topology. The widgets allow to specify:\n", " * `netParam`: the linking probability with which each agent has a communication edge with another agent\n", " * `barabasi-albert`: Random network with Barabasi-Albert topology. The widgets allow to specify:\n", " * `netParam`: the number of edges every new node (during creation) is linked to\n", " * `dynamic`: Simulation of moving particles in which agent move in a periodic space and interact locally. The widgets allow to specify:\n", " * `netParam`: range of interaction of the agent (radius)\n", " * `motionCorrelatedness`: motion correlatedness of the agent. Correlatedness = 1 correponds to straight movement, Correlatedness = 0 corresponds to uncorrelated random walk.\n", " * `particleSpeed`: speed of the agent, i.e. displacement in one timestep\n", " * `showTrace`: flag to show the motion trajectory of the agents\n", " * `showInteractions`: flag to show interactions; agents within communication range are liked through a line\n", "* `timestepSize`: length of one timestep; the maximum timestep is determined by the rates\n", "\n", "(Note that in addition to the ones specified for the `SSA` simulations, here the widgets to show motion trajectory and to show interactions do not trigger the recomputation of the simulations).\n", "Static networks are available only for models that do not contain the empty-set or constant reactant, whereas moving particles have no constraints; see [Supported features](#features)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Advanced options for `integrate()` and `noiseCorrelations()`\n", "* `initialState`: specify the initial conditions for all reactants in the systems. Note that if `conserved = True` occurs as a keyword argument in `integrate()` and `noiseCorrelations()` then initial conditions cannot be set independently. \n", "* `maxTime`: specify the length of the simulation.\n", "* `plotProportions` (only in `integrate()`): choose if absolute numbers or proportions are displayed (check box). Proportions are denoted $\\Phi$ with the reactant as index. If `True`, the absolute number is divided by the total number of reactants (system size) at $t=0$. If the system is conserved, the sum of all concentrations sum to 1 for all $t$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "#### Keyword arguments\n", "\n", "\n", "\n", "\n", "##### Keyword arguments for `stream()` and `vector()`\n", "\n", "- valid for both methods:\n", "\n", " - `showFixedPoints = False`: plot fixed points\n", " \n", " - `showNoise = False`: plot noise around fixed points\n", "\n", " - `runs = 1`: number of simulation runs to be executed, must be strictly positive\n", "\n", " - `aggregateResults = True`: aggregate the results from multiple runs\n", " \n", " - `fontsize = None`: specify fontsize for axis-labels, accepts integers\n", " \n", " - `xlab = 'stateVariable1'`: specify label on x-axis\n", " \n", " - `ylab = 'stateVariable2'`: specify label on y-axis\n", " \n", " - `choose_xrange = None`: specify range plotted on x-axis; usage: `choose_xrange=[xmin, xmax]`, `xmin` and `xmax` are of type: float, if not given uses data values to set axis limits\n", " \n", " - `choose_yrange = None`: specify range plotted on y-axis; usage: `choose_xrange=[ymin, ymax]`, `ymin` and `ymax` are of type: float, if not given uses data values to set axis limits\n", " \n", " - `silent = False`: switch off widgets and plot, important for use with multi-controllers\n", " \n", "- valid only for `vector()`:\n", "\n", " - `zlab = 'stateVariable3'`: specify label on z-axis, 3D plots only (**N.B.** 3D stream plots not yet implemented)\n", " \n", "##### Keyword arguments for `bifurcation()`\n", "\n", " - `initialState = {}`: initial proportions of the reactants (type: float in range [0,1]), will be used ONLY if calculation of stationary states fails, can also be set via `initWidgets` argument\n", " \n", " - `initBifParam = 2.0`: initial value of bifurcation parameter, can also be set via `initWidgets` argument\n", " \n", " - `conserved = False`: specify if a system is conserved to make proportions of state variables (at time t=0) sum up to 1, or not; type: boolean (takes True or False)\n", " \n", " - `contMaxNumPoints = 100`: choose maximum number of continuation points, accepts integers\n", " \n", " - `fontsize = None`: specify fontsize for axis-labels, accepts integer, if not given fontsize is automatically derived from length of axis label\n", " \n", " - `xlab = 'bifurcationParameter'`: specify label on x-axis (type: str), if not given uses symbol for `bifurcationParameter` in arguments as default label\n", " \n", " - `ylab = r'\\Phi_{stateVariable1}'`: specify label on y-axis (type: str), if not given uses expression given as `stateVariable1` in arguments as index, if `stateVariable1` is a sum/difference of the form $Reactant_1 +/- Reactant_2$ the default label is $\\Phi_{Reactant1} +/- \\Phi_{Reactant2}$, i.e. expressions will be LaTeX-rendered, hence the `r` before the string which converts it into raw string literals (remeber we use $\\Phi$ to denote concentrations)\n", " \n", " - `choose_xrange = None`: specify range plotted on x-axis; usage: `choose_xrange=[xmin, xmax]`, xmin and xmax are of type: float, if not given uses data values to set axis limits\n", " \n", " - `choose_yrange = None`: specify range plotted on y-axis; usage: `choose_xrange=[ymin, ymax]`, ymin and ymax are of type: float, if not given uses data values to set axis limits\n", " \n", " - `silent = False`: switch on/off widgets and plot (type: boolean), accepts `True` or `False`; important for use with multi-controllers\n", "\n", "##### Keyword arguments for `SSA()` and `multiagent()`\n", "\n", "- valid for both methods:\n", "\n", " - `params = None`: list of parameters defined as pairs ('parameter name', value), see [Partial controllers](#partialcont). Rates defaults to `mumot.MuMoTdefault._initialRateValue`, system size defaults to `mumot.MuMoTdefault._systemSize`\n", "\n", " - `initialState = None`: initial proportions of reactants, dictionary with reactants as keys and floats in range [0,1] as values. Defaults to all reactants being in alphabetically first state. See [bookmarks](#bookmarking).\n", "\n", " - `maxTime = mumot.MuMoTdefault._maxTime`: simulation time, must be strictly positive\n", " \n", " - `randomSeed = random`: random seed, must be strictly positive in range [0, `mumot.MAX_RANDOM_SEED`]\n", " \n", " - `plotProportions = False`: plot proportions otherwise plot full populations, defaults to `False`\n", " \n", " - `realtimePlot = False`: plot results in realtime\n", " \n", " - `visualisationType = 'evo'`: type of visualisation (``'evo'``,``'final'`` or ``'barplot'``)\n", " \n", " - `final_x`: which reactant is shown on x-axis when visualisation type is ``'final'``. Defaults to the alphabetically first reactant\n", " \n", " - `final_y`: which reactant is shown on y-axis when visualisation type is ``'final'``. Defaults to the alphabetically second reactant\n", "\n", " - `runs = 1`: number of simulation runs to be executed, must be strictly positive\n", "\n", " - `aggregateResults = True`: aggregate the results from multiple runs, accepts True or False\n", " \n", " - `legend_fontsize = 14`: specify fontsize of legend (type: integer)\n", "\n", " - `legend_loc = 'upper left'`: specify legend location: combinations like 'upper left', 'lower right', or 'center center' are allowed (9 options in total)\n", "\n", " - `fontsize = None`: specify fontsize for axis-labels, accepts integer, if not given fontsize is automatically derived from length of axis label\n", "\n", " - `xlab = 'time t'`: specify label on x-axis (type: str), if given replaces default\n", "\n", " - `ylab = 'noise correlations'`: specify label on y-axis (type: str), if given replaces default\n", "\n", " - `choose_xrange = None`: specify range plotted on x-axis; usage: `choose_xrange=[xmin, xmax]`, xmin and xmax are of type: float, if not given uses data values to set axis limits\n", "\n", " - `choose_yrange = None`: specify range plotted on y-axis; usage: `choose_xrange=[ymin, ymax]`, ymin and ymax are of type: float, if not given uses data values to set axis limits\n", "\n", " - `silent = False`: switch off widgets and plot, important for use with multicontrollers\n", "\n", "\n", "- valid only for `multiagent()`:\n", "\n", " - `netType = 'full'`: type of network. Available types are: `full`, `erdos-renyi`, `barabasi-albert`, `dynamic` (see more information about their meaning above)\n", " - `netParam`: this parameters has a specific meaning depending on the `netType`, see details above\n", " - `motionCorrelatedness = 0.5`: (only for `netType`== `dynamic`) motion correlatedness of the agent. Correlatedness = 1 correponds to straight movement, Correlatedness = 0 corresponds to uncorrelated random walk.\n", " - `particleSpeed = 0.01`: (only for `netType`== `dynamic`) speed of the agent, i.e. displacement in one timestep\n", " - `showTrace = False`: (only for `netType`== `dynamic`) flag to show the motion trajectory of the agents\n", " - `showInteractions = False`: (only for `netType`== `dynamic`) flag to show interactions; agents within communication range are liked through a line\n", " - `timestepSize`: length of one timestep; the maximum timestep is determined by the rates\n", "\n", "##### Keyword arguments for `integrate()` and `noiseCorrelations()`\n", "\n", "- valid for both methods:\n", "\n", " - `maxTime = 3.0`: simulation time for temporal evolution of reactants in `integrate()` or correlations of noise in `noiseCorrelations()`, type: float (larger than 0)\n", "\n", " - `tstep = 0.01`: time step used for numerical integration of reactants in `integrate()` or correlations of noise in `noiseCorrelations()`, type: float (larger than 0)\n", "\n", " - `initialState = {}`: initial proportions of the reactants (type: float in range [0,1]), can also be set via `initWidgets` argument\n", "\n", " - `conserved = False`: specify if a system is conserved to make proportions of state variables (at time t=0) sum up to 1, or not; type: boolean (takes True or False)\n", "\n", " - `legend_fontsize = 14`: specify fontsize of legend (type: integer)\n", "\n", " - `legend_loc = 'upper left'`: specify legend location: combinations like 'upper left', 'lower right', or 'center center' are allowed (9 options in total)\n", "\n", " - `fontsize = None`: specify fontsize for axis-labels, accepts integer, if not given fontsize is automatically derived from length of axis label\n", "\n", " - `xlab = 'time t'`: specify label on x-axis (type: str), if given replaces default\n", "\n", " - `ylab = 'noise correlations'`: specify label on y-axis (type: str), if given replaces default\n", "\n", " - `choose_xrange = None`: specify range plotted on x-axis; usage: `choose_xrange=[xmin, xmax]`, xmin and xmax are of type: float, if not given uses data values to set axis limits\n", "\n", " - `choose_yrange = None`: specify range plotted on y-axis; usage: `choose_xrange=[ymin, ymax]`, ymin and ymax are of type: float, if not given uses data values to set axis limits\n", "\n", " - `silent = False`: switch on/off widgets and plot (type: boolean), accepts True or False; important for use with multi-controllers\n", "\n", "\n", "- valid only for `integrate()`:\n", "\n", " - `plotProportions = False`: flag to plot proportions or full populations (type: boolean)\n", " \n", " \n", "- valid only for `noiseCorrelations()`:\n", "\n", " - `maxTimeDS = 50`: simulation time for ODE system (type: float larger than 0); the final state of the ODE-system after `maxTimeDS` is the state around which noise correlations are calculated and plotted (MuMoT checks if this is a stationary state)\n", " \n", " - `tstepDS = 0.01`: time step of numerical integration of ODE system (type: float)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Advanced users" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Seeing verbose error messages\n", "\n", "By default, error messages are brief to avoid distracting the majority of users with unnecessary details. However, if you are an advanced user or developer, you may need to see these details. To make error messages more verbose and enable *traceback* run the following command..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mumot.setVerboseExceptions(True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Downloading results\n", "\n", "Some views allow access to their raw results, such as time series from stochastic simulations. To access these and analyse in an external application, click on the *save file* button next to the bookmark button; a `Download link` will appear which can be clicked to open the raw results in a new browser window. Note that not all views allow results download at present (see [Supported features](#downloadI)).\n", "\n", "If using a view independently of a controller, if it implements the `downloadResults()` method this also produces a download link." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ssa2 = model1.SSA()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ssa2._view.downloadResults()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Partial controllers\n", "\n", "If you are creating a notebook to demonstrate a concept, it may be easiest to fix some of the parameters to ensure that the user can only vary the parameters you want them to. Partial controllers can be created by passing the standard commands a list `params` of parameter name, parameter value pairs. For example..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "partial1 = model4.stream('A','B', params = [('\\\\mu',1.5),('s',1)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Initialising widgets\n", "\n", "Widget initial values, as well as ranges and step sizes where applicable, can all be specified. For example..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "init1 = model4.SSA(initWidgets = {'\\Delta':[0, -2, 2, 0.2], 'realtimePlot':True, 'visualisationType':'barplot'})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Multicontrollers\n", "\n", "Since MuMoT provides you with multiple views on the same model, it can be useful to have a single controller updating the parameters of multiple views simultaneously; for example, we may wish to compare stochastic simulation results against the infinite-population ODE solutions. MuMoT can do this, but since this is advanced functionality the syntax is a little more complicated.\n", "\n", "\n", "We begin by making 'silent' versions of the controllers, using the keyword `silent = True`\n", "\n", "Then we combine these by building a new instance of the `MuMoTmultiController` class and passing in the list of constituent controllers we want to use. If the views should all be plotted on the same axes we use the keyword `shareAxes = True` (**N.B.** the legend to be used will be taken from the last controller in the list). Widgets can be initialised using the `initWidgets` keyword (see the section [Keyword arguments for controllers](#kwargscontrollers)); in general all appropriate keywords can be sent to the multicontroller." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "silentIntegrate1 = model4.integrate(silent = True)\n", "silentSSA1 = model4.SSA(silent = True)\n", "mc1 = mumot.MuMoTmultiController([silentIntegrate1, silentSSA1], shareAxes = True,\n", " initWidgets={'initialState':{'U': [1,0,1,0.01],'B': [0,0,1,0.01],\n", " 'A': [0,0,1,0.01]}, 'runs':[5,1,20,1], 'aggregateResults':True})" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mc1.showLogs(tail = True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is possible to construct nested multicontrollers, by making the component multicontrollers silent. For example, making a silent version of the previous multicontroller, which can then be combined with another (multi)controller..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "silentIntegrate2 = model4.integrate(silent = True)\n", "silentSSA2 = model4.SSA(silent = True)\n", "mc2 = mumot.MuMoTmultiController([silentIntegrate2, silentSSA2], shareAxes = True, silent = True, \n", " initWidgets={'initialState':{'U': [1,0,1,0.01],'B': [0,0,1,0.01],\n", " 'A': [0,0,1,0.01]}, 'runs':[5,1,20,1], 'aggregateResults':True})\n", "\n", "mc3 = mumot.MuMoTmultiController([mc2, model4.stream('A', 'B', silent = True)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Of course, all multicontrollers can be [partial](#partialcont) and [bookmarked](#bookmarking)..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "silentIntegrate3 = model4.integrate(silent = True)\n", "silentSSA3 = model4.SSA(silent = True)\n", "silentStream1 = model4.stream('A', 'B', silent = True)\n", "mc4 = mumot.MuMoTmultiController([silentIntegrate3, silentSSA3], shareAxes = True, silent = True, \n", " initWidgets={'initialState':{'U': [1,0,1,0.01],'B': [0,0,1,0.01],\n", " 'A': [0,0,1,0.01]}, 'runs':[5,1,20,1], 'aggregateResults':True})\n", "\n", "mc5 = mumot.MuMoTmultiController([mc4, silentStream1], params = [('\\mu',1.5), ('s',1)])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mc5.showLogs(tail = True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Checking the stoichiometry of a model's reactions\n", "\n", "\n", "In MuMoT the stoichiometry of a *reaction* characterising the relationship between *reactants* and *products* can be accessed via `getStoichiometry()` and displayed using `showStoichiometry()`. The method `getStoichiometry()` allows you to access a model-specific dictionary with its keys being the *reaction numers* (*reaction number* represents another dictionary with reaction rate, number of particles before and after the reaction, and substitutions of state variables if applicable). Here is an example using our `model4`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model4.getStoichiometry()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Whereas this output is useful if the items of the dictionary need to be accessed,\n", "`showStoichiometry()` produces nicely rendered output (whilst containing the same information):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model4.showStoichiometry()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Viewing the Master and Fokker-Planck equations, and calculating & displaying noise\n", "\n", "\n", "Just to make sure you really want to look at these, and because the size of the equations can slow down your browser, some of the following commands are commented. Uncomment then run if you want to see the gory details..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To show this functionality we make use of a model describing the production and decay of protein $P$ and its dimerization into $P_2$ (more conveniently called $X$ and $Y$ below), see Hayot & Jayaprakash ([2004](#references)) for details. This model is defined as:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%model\n", "$\n", "\\emptyset -> X : k_3\n", "X -> \\emptyset : k\n", "X + X -> Y + \\emptyset : k_1\n", "Y + \\emptyset -> X + X : k_2\n", "$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you already know from the examples above, first we must parse the model ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model7 = mumot.parseModel(In[-2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that MuMoT has parsed the model, we can continue deriving the *chemical Master equation*. \n", "\n", "Master equations describe the time evolution of a system that, in general, can be in one of many different possible states at a given point in time with a certain probability. Over time, the system makes transitions so that the overall state changes continuously. Transitions are governed by transition probabilities that are usually summarised in a transition matrix. The Master equation is a set of first-order differential equations that describe the evolution of the state-dependent probabilities driven by the transition matrix.\n", "\n", "The chemical Master equation ([van Kampen, 2007](#references)) is a Master equation for the probabilities that the system has any given composition of rectants as a function of time. In MuMoT this equation can be derived using:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model7.showMasterEquation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, $P(X,Y,t)$ describes the probability of observing the systems's composition of reactants $X$ and $Y$ at time $t$. In MuMoT, the system size that, at this point, is not specified and what is used as a formal parameter for the system size expansion described below is denoted $\\overline{V}$. Please pay attention to this notation if defining your model using the letter $V$ (but don't worry, MuMoT will not confuse both variables).\n", "\n", "Here *step operators* $E_{op}$ are used to express the chemical Master equation in a compact form. MuMoT implements the derivation of the Master equation as shown by van Kampen ([2007](#references)).\n", "\n", "Unfortunately, there is only a limited number of simple cases that allow exact solution of the Master equation. Hence, analytical approximations may help, such as the *system size expansion* (also called *van Kampen expansion*), which is implemented in MuMoT by `showVanKampenExpansion()`. \n", "\n", "The idea of this method is to expand the Master equation in powers of a small parameter, neglecting all except the leading term, which is\n", "
\n", "$$\\sim \\frac{1}{\\sqrt{system\\, size}} \\,.$$\n", "
\n", "\n", "Let's do the expansion with MuMoT and display it for the current model ($\\overline{V}$ = system size):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model7.showVanKampenExpansion()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This equation describes an expansion in the parameter $\\overline{V}$, which represents the system size. We repeat that at this point this is not specified. However, often $\\overline{V}$ is a volume or the total number of particles in a real system. For more information on the system size expansion refer to van Kampen ([2007](#references)).\n", "\n", "Collecting all terms $\\propto \\sqrt{\\overline{V}}$ gives the ODEs which, as described at the very beginning of this notebook, can be displayed using" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model7.showODEs(method='vanKampen')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, the $\\Phi_{X,Y}$ once more describe the concentrations of $X$ and $Y$ in the mixture (i.e. $\\Phi_{X}=X/\\overline{V}$, for example).\n", "\n", "Next, if we collect all terms $\\propto \\overline{V}$ we obtain the *Fokker-Planck equation* in the *linear noise approximation*:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model7.showFokkerPlanckEquation()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This equation describes the probability of fluctuations $\\eta_{X,Y}$ at time $t$. In MuMoT, fluctuations are automatically denoted $\\eta$ with the index denoting the reactant. These quantities describe the fluctuations around the *mean field* value of the concentrations $\\Phi_{X,Y}$. \n", "\n", "Having this equation allows us derive the equations of motion for the fluctuations. This can be done as follows:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model7.showNoiseEquations()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, we have first-order moments, $\\langle \\eta_{X}\\rangle$ and $\\langle \\eta_{Y}\\rangle$, and second-order moments, $\\langle \\eta_{R_1}\\eta_{R_2}\\rangle$, where $R_{1,2}=X,Y$.\n", "\n", "We can analytically derive the stationary solutions of these equations and display them using" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "model7.showNoiseSolutions()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another possibility to study noise in the system is by calculating the time-dependent noise correlations according to" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nc1 = model7.noiseCorrelations(initWidgets={'maxTime':[20,5,50,1]})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These correlation functions relate to fluctuations around a fixed point. To see the corresponding fixed point we can simply show the logs:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nc1.showLogs()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Please keep this stationary state of the system in mind for a little moment.\n", "\n", "We can also numerically integrate our dynamical system: " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "int2 = model7.integrate(initWidgets={'maxTime':[30,5,50,1]}, legend_loc='center right')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now click on the *Advanced options* tab above the figure and tick the *Plot population proportions* tab. We then display the last point on the curve by showing the logs with the following command:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "int2.showLogs()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, you can see the last point on the curve before selecting to plot proportions and after choosing this option. Remember, the proportions are calculated by dividing the total numbers by the system size at $t=0$ (which, if you haven't changed anything, should equal 10 in the plot above -- see the *System size* slider above the figure). \n", "\n", "Note that the values on the sliders are the same as in the plot of the noise correlations above. If you can recall the fixed point around which the correlation functions were plotted (if not, it was ($\\Phi_X=1$, $\\Phi_Y=1$)) and compare this with the last point on the curve above (see `int2.showLogs` above), we can safely say that both are in agreement. So, it works! The noise correlations have been plotted around the stationary state of the dynamical system. Using `integrate` we obtained a graphical confirmation of this." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Accessing equation objects\n", "\n", "Expert users can gain access to symbolic equation objects representing all the derived equations within MuMoT (see [Command reference](#comreference). These are [Sympy](https://www.sympy.org/en/index.html) objects, suitable for further analysis and manipulation. A list of equations is provided; if a system-size substitution has been performed then the equation for this is appended to the list; for example" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "eqlist = model2.getVanKampenExpansion()\n", "eqlist[2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is beyond the scope of this user manual at present to provide detailed guidance on using these objects." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### References\n", "\n", "* Gillespie, D. T. (1976). A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. *Journal of Computational Physics* **22** (4): 403–434\n", "* Hayot, F. & C. Jayaprakash (2004) \n", "[The linear noise approximation for molecular fluctuations within cells](http://iopscience.iop.org/article/10.1088/1478-3967/1/4/002/pdf). *Physical Biology* **1**, 205-210\n", "* Marshall, J. A. R., Reina, A. Bose, T. (2019) Multiscale Modelling Tool: Mathematical modelling of collective behaviour without the maths. *In submission*\n", "* Murray, J. D. (2002) *Mathematical Biology I. An Introduction (Second Edition)*. Springer\n", "* Pais, D., Hogan, P.M., Schlegel, T., Franks, N.R., Leonard, N.E. & Marshall, J.A.R. (2013) [A mechanism for value-sensitive decision-making](http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0073216). *PLoS one* **8**(9), e73216\n", "* Reina A., Valentini G., Fernández-Oto A., Dorigo M. & Trianni V. (2015) [A design pattern for decentralised decision making](http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0140950). *PLoS one* **10**(10), e0140950\n", "* Reina A., Marshall J.A.R., Trianni V., Bose T. (2017) [Model of the best-of-N nest-site selection process in honeybees](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.95.052411) *Physical Review E* **95**: 052411\n", "* Seeley, T.D, Visscher, P.K. Schlegel, T., Hogan, P.M., Franks, N.R. & Marshall, J.A.R. (2012) [Stop signals provide cross inhibition in collective decision-making by honeybee swarms](http://www.sciencemag.org/content/335/6064/108.full.pdf). *Science* **335**, 108-111\n", "* van Kampen, N. (2007) *Stochastic Processes in Physics and Chemistry (Third Edition)*. North-Holland" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.4" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": { "height": "calc(100% - 180px)", "left": "10px", "top": "150px", "width": "253px" }, "toc_section_display": true, "toc_window_display": true }, "widgets": { "state": { "0f1da6289ccc47ec96c8c8b477f3e37e": { "views": [ { "cell_index": 60 } ] }, "12dff49951bf4f4aa6afd2437d2733ea": { "views": [ { "cell_index": 38 } ] }, "143c3aa5b3a94602abcb99f5d3682b42": { "views": [ { "cell_index": 23 } ] }, "15f8da4f248e41fdbb86e469bbf0a423": { "views": [ { "cell_index": 51 } ] }, "22efce659af24047800463f2c6d460ee": { "views": [ { "cell_index": 23 } ] }, "2305e4265a234613af7ec1488c3a84a0": { "views": [ { "cell_index": 51 } ] }, "25eaa6a79b394d96bf6b016a8f5789a0": { "views": [ { "cell_index": 33 } ] }, "25f7cb6ad4594dbb9c3dfaeeeebd2baa": { "views": [ { "cell_index": 33 } ] }, "2a889abb13a441a6b0b1bb07b06aae14": { "views": [ { "cell_index": 60 } ] }, "2db99fe905e1430a9e0fc3f80e61824d": { "views": [ { "cell_index": 33 } ] }, "31ffe852616040109a3b4f88a6815b8a": { "views": [ { "cell_index": 36 } ] }, "3e3b2977642f48e99ceeef154551af71": { "views": [ { "cell_index": 60 } ] }, "4ba771e07fdb44c49a07b94aafee0c1d": { "views": [ { "cell_index": 38 } ] }, "4d722defa6f54a03956ef26e5bb24d75": { "views": [ { "cell_index": 38 } ] }, "4e4fb5a6665b4e89acf172973b5f9b27": { "views": [ { "cell_index": 23 } ] }, "50ab0a9b7d4d47d2bef8033efa6b6661": { "views": [ { "cell_index": 38 } ] }, "53a25e884ed74eaab9813cb4e3413c2f": { "views": [ { "cell_index": 38 } ] }, "56f18bd9baa1474fac1d692af347424f": { "views": [ { "cell_index": 51 } ] }, "5f4fd39c4f5047db8441d0b30df86a38": { "views": [ { "cell_index": 38 } ] }, "68117c20182148de957ae1a311729bc2": { "views": [ { "cell_index": 23 } ] }, "6b9e63ee22924ae3b6731d530ba6b549": { "views": [ { "cell_index": 36 } ] }, "7050fe3f432e4180b329a3c92ed09551": { "views": [ { "cell_index": 51 } ] }, "7713dfc8cf164108b44b49d477f3e154": { "views": [ { "cell_index": 51 } ] }, "79f7fdc5299e4b1bb75d32bb0e87d7a3": { "views": [ { "cell_index": 23 } ] }, "7b957930c9924cb3b1cca132418e2f00": { "views": [ { "cell_index": 51 } ] }, "8b82dea0b0744375b10edd290ea509d3": { "views": [ { "cell_index": 36 } ] }, "8d07f6b41a444b49b12d7f8e5918ca6a": { "views": [ { "cell_index": 23 } ] }, "8ef0ad62478b43129d8cd652add052c7": { "views": [ { "cell_index": 23 } ] }, "9560997a50274839a7117e95b073afba": { "views": [ { "cell_index": 36 } ] }, "96b329f0ab1d47338ac6d63da570c275": { "views": [ { "cell_index": 36 } ] }, "9adb34f1a17041719156750646d42cfb": { "views": [ { "cell_index": 51 } ] }, "9e7afa9d8d1b478a83a0955428961503": { "views": [ { "cell_index": 23 } ] }, "a408184986ae423dbde56361209f036a": { "views": [ { "cell_index": 33 } ] }, "a7ca5685289644b8aa449ce9377e6a5c": { "views": [ { "cell_index": 23 } ] }, "a7dbac4a99a04c4e931b4cc6c9beda28": { "views": [ { "cell_index": 51 } ] }, "a874f3c5de9548deab1f15d01e0a35aa": { "views": [ { "cell_index": 60 } ] }, "ab5a6f55976e4b99847157fde376160d": { "views": [ { "cell_index": 38 } ] }, "ad3aa76f7e874016acdd53204362e367": { "views": [ { "cell_index": 51 } ] }, "b64577ba4bc14887b39249c7599693e2": { "views": [ { "cell_index": 33 } ] }, "b6818d31691542d2b4fc017e2aad1400": { "views": [ { "cell_index": 23 } ] }, "c6e9aa6f9e6b4bbf9e994e85216e6cbf": { "views": [ { "cell_index": 33 } ] }, "caac276fb0844d3b939ae00b7ae0ab95": { "views": [ { "cell_index": 23 } ] }, "d8aaa2637f8b4586b243e7d1e265614a": { "views": [ { "cell_index": 51 } ] }, "e77472d2a0524ea0b75dcee2471685e4": { "views": [ { "cell_index": 60 } ] }, "f2876ad0e37e4b748ff80f9cc3791542": { "views": [ { "cell_index": 23 } ] }, "f73cd695eb284175a7c65821b942925c": { "views": [ { "cell_index": 23 } ] }, "f8319629082441e98a925c1d162f75c5": { "views": [ { "cell_index": 60 } ] }, "f8f9d7c2ae4d4262b564c5f434f7f2d0": { "views": [ { "cell_index": 33 } ] }, "fd6fcf6218cc4d3bb180fbb1264c115d": { "views": [ { "cell_index": 36 } ] } }, "version": "1.2.0" } }, "nbformat": 4, "nbformat_minor": 2 }