Skip to content

Tutorial 1

Tim Vaughan edited this page Apr 4, 2014 · 34 revisions

Stochastic Logistic Model Simulation

In this first tutorial, we will demonstrate a very basic use of MASTER by using it to generate a simulated population history under a stochastic version of the logistic model. We will also demonstrate how to visualize the result using the statistical analysis package R.

Specifically, the stochastic model we will simulate under involves a population of individuals of a single type X evolving under the continuous time birth-death process described by the following pair of reactions:

The first may describe asexual reproduction of an individual and occurs at an average rate while the second describes a non-linear death process as may be expected to occur as a result of resource competition and which occurs at the average rate . The model is so-named because in the large population size limit, the expected population size n(t) obeys the logistic growth equation:

We thus expect the population size of the stochastic model to stabilize near the deterministic carrying capacity for any initial population size n(t) >= 1.

Input file specification

We will now piece together a MASTER input file describing the model and other information needed to generate the simulated population size history. Firstly, we open a blank text file in a text editor. We then enter the following XML element which serves to (a) notify the program that this is a valid input file and (b) specify the locations of the BEAST plugins constituting MASTER. (This element must exist in every MASTER input file.)

<beast version='2.0' namespace='master:master.model:master.steppers:master.conditions:master.postprocessors:master.outputs'>
</beast>

Within this element, we then add the <run> element which specifies high-level details about the type of calculation to be run. In our case, we desire a single simulated trajectory with a fixed length of 20 days. The XML thus becomes:

<beast version='2.0' namespace='master:master.model:master.steppers:master.conditions:master.postprocessors:master.outputs'>
    <run spec='Trajectory'
         simulationTime='20'>

    </run>
</beast>

Model

We now come to the specification of the model itself. This is achieved by placing a model element within <run>. In the case of the stochastic logistic model described above, the <model> element takes the form

<model spec='Model' id='model'>
    <population spec='Population' populationName="X" id='X'/>
    <reaction spec='Reaction' reactionName="Birth" rate="1.0">
        X -> 2X
    </reaction>
    <reaction spec='Reaction' reactionName="Death" rate="0.01">
        2X -> X
    </reaction>
</model>

This model contains two distinct element types. Firstly, the population element is used to indicate that the model includes an unstructured population named "X". (The value of the id attribute is distinct from this name and is used to refer to this particular population element from elsewhere in the XML file.) Secondly, the reaction elements specify the individual birth and death processes in the model in a very similar form to that given in the reactions at the start of this page. Each of these processes/reactions are given an appropriate name and rate via the reactionName and rate attributes. In this case we have chosen rates and , corresponding to a carrying capacity of 100 in the deterministic limit.

Initial conditions

An equally important ingredient in specifying the details of the simulation is the initial condition, specified via the initialState element. For simple population-level stochastic simulations such as this, initial state specification equates simply to specifying the initial size of the constituent populations. These sizes are given by adding populationSize elements as necessary. In the case of the stochastic logistic model we only need one of these elements:

<initialState spec='InitState'>
    <populationSize spec='PopulationSize' population="@X" size='1'/>
</initialState>

Note that a population attribute is used to specify the particular population that the size refers to. In this case, we have used a BEAST reference to the population element given earlier in the XML. Such references are made up of an "@" followed by the BEAST ID of the element we're referring to, which in this case is "X".

Note that populations not explicitly set to an initial value will be implicitly assumed to have an initial value of zero.

Output specification

The final component of our input file is involved in determining what MASTER does with the simulated history it generates. This is achieved through the use of the following output element:

<output spec='JsonOutput' fileName='StochasticLogistic_out.json'/>

which tells MASTER that the output should be written in JSON format to the file StochasticLogistic_out.json.

Full input file

Bringing all of these elements together, we have the final MASTER input file specifying the required calculation:

<beast version='2.0' namespace='master:master.model:master.steppers:master.conditions:master.postprocessors:master.outputs'>
    <run spec='Trajectory'
         simulationTime='20'
         verbosity='1'>

        <model spec='Model' id='model'>
            <population spec='Population' populationName="X" id='X'/>
            <reaction spec='Reaction' reactionName="Birth" rate="1.0">
                X -> 2X
            </reaction>
            <reaction spec='Reaction' reactionName="Death" rate="0.01">
                2X -> X
            </reaction>
        </model>

        <initialState spec='InitState'>
            <populationSize spec='PopulationSize' population="@X" size='1'/>
        </initialState>
        
        <output spec='JsonOutput' fileName='StochasticLogistic_out.json'/>
    </run>
</beast>

Running MASTER

The XML above should be saved in a file named something sensible such as StochasticLogistic.xml. The specifics of how MASTER is executed will depend on your operating system and are described on the main page. Once this process is complete, however, the result will be the same - a file named StochasticLogistic_out.json containing the results of the simulation will appear ready for additional processing or visualization.

Output visualization using R

Loading the simulated history into R requires use of the rjson package which is available from CRAN. This can be installed directly from within R by issuing the following command at the prompt:

> install.packages('rjson')

Once this command is entered, a dialog box may appear asking the user to select from one of the available mirrors. After doing this, rjson will be automatically downloaded and installed. It must then be loaded using the following command:

> library(rjson)

Once rjson is installed and loaded, we need to ensure that R's working directory matches the location of the MASTER output file. This is done using the following setwd command, with '/path/to/input/' replaced by the absolute path of the directory containing the output file.

> setwd('/path/to/output/')

To actually load the file, we use the fromJSON function. Note that the file name must be given as a named argument. (If this name is omitted, the function will attempt to parse the file name itself as a JSON string.) The following command loads the simulation output and stores it in in the variable df.

> df <- fromJSON(file='StochasticLogistic_out.json')

Once loaded, we can interrogate the data variable using the names function:

> names(df)
[1] "t"   "sim" "X"  

Here we see that df contains three distinct elements. The element df$t contains an array of the times at which the simulated population size has been recorded. The population sizes themselves are contained in the array df$X. The remaining element df$sim contains a variety of information about the calculation itself and can be useful for keeping track of the conditions under which the output was generated.

The simulated population history can now be visualized using the following plot command, which should open up a new plot window illustrating the results. (If you have problems getting the plot window to appear, you should consult the R documentation specific to your system.)

> plot(df$t, df$X, 's', col='blue', xlab='Time', ylab='Population size')
> lines(c(0,30), c(100,100), lty=2, col='red', lwd=2)

The result of this command is shown below. The carrying capacity from the deterministic approximation is also shown for comparison.