Tutorial 1
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.
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>
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.
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.
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
.
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>
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.
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.