MDSL Documentation

Membrane tags are defined in square brackets before the membrane type but after the number.

There are two different ways to refer to a membrane within an MDSL file: by its type name and by one of its tags. The purposes of these two ways are:

• Types allow different reactions to happen in different membrane types
• Tags allow species to be initialised differently in different membranes of the same type

There are two different ways in which tags can be used:

1. To initialise different instances of the same membrane type differently
2. To initialise membranes of different types the same

Initialising different instances of the same membrane type differently

                                initial tree { simulation { 19 [not-infected] cell } { 1 [infected] cell } }
                                species pathogen contained infected = 100 units
                            

Initialising membranes of different types the same

                                initial tree { simulation {3 [t] A} {4 B} {2 [t] C} }
                                species x contained t = 30 units
                                species y contained t = 12 units
                                species z contained t = 136 units
                            

                            species r on cell = 10 receptors.cell^-1
                        

Species definitions set the initial number of each species present at the start of the simulation.

The units must be specified in each species definition, but these units are not used or checked during the simulation.

If a species definition references a membrane that there are multiple copies of in the initial tree, then that number of the species is added to each membrane. For example, the following code creates three copies of membrane A, each with five copies of species x inside. Thus there will be 15 species x in the whole simulation.

                            initial tree {simulation {3 A}}
                            species x contained A = 5 units  # 15 copies of species x created
                        

NB: for species definitions that are specified with an around location, the around will resolve to the parent membrane before the decision of how many species to add. So the following code will create only three of species x.

                            initial tree {simulation {3 A}}
                            species x around A = 5 units  # 3 copies of species x created
                        

                            a_r_complex on cell modifier k_signal => a_r_complex on cell and 2 b contained cell
                        

The left hand side (reactants) and right hand side (products) of reactions are species in locations.

                            a contained A modifier 3 * (x contained A + y contained A) => b contained A
                        

Reaction modifiers are equations, which can contain located species as terms in the equation.

When used as terms in a modifier equation, the species numbers will not be modified when the reaction fires, but they will influence the rate of the reaction. This can be used to create rate equations that would not be possible by just including these species in the stoichiometry of the reaction, for example the reaction above, with modifier 3 * (x contained A + y contained A).

It can be cumbersome to write the locations for each species all the time, particularly when species are used in reaction modifiers.

If multiple species are referenced within the same location, then a location context block can be used to simplify the MDSL.

For example, the following two snippets of code define the same model:

                            species a contained M = 1 units
                            species b contained M = 2 units
                            species c on M = 3 units

                            parameter k = 4 units

                            a contained M binds b contained M modifier k => b contained M
                            c on M modifier (a contained M + b contained M) / b contained M => d under M
                        

                            contained M {
                              species a = 1 units
                              species b = 2 units
                              species c on M = 3 units

                              parameter k = 4 units

                              a binds b modifier k => b
                              c on M modifier (a + b) / b => d under M
                            }
                        

Location context blocks automatically insert a given location into any species names that are used without a location inside the block. Other locations can be used inside the block (for example on M and under M above), as long as they are specified explicitly.

To aid in checking that parameter equations have been written correctly, the simulator outputs a file named output_Parameters.txt, containing all the parameter definitions, along with the numbers that they evaluate to. An example of this file is shown below:

                            Parameters equations defined in the model file:
                            parameter k_on = (0.2) ml.ng^-1.cell.receptors^-1.hour^-1
                            parameter k_d = (0.005) ng.ml^-1
                            parameter k_off = ((k_on) / (k_d)) cell.receptors^-1.hour^-1

                            Parameter values computed:
                            parameter k_on = 0.2 ml.ng^-1.cell.receptors^-1.hour^-1
                            parameter k_d = 0.005 ng.ml^-1
                            parameter k_off = 40.0 cell.receptors^-1.hour^-1
                        

The first part of this file shows the parameter equations with all their implicit brackets added. This makes clear the order of precedence of the operators in the equations.

The second part of this file shows the actual numbers that the equations evaluate to.

The reaction modifier is the equation that comes immediately after the modifier keyword.

This can be a simple rate constant: a number or a parameter (in the above example, it is the parameter k). Or it can be a complex equation, potentially involving multiple parameters and species concentrations.

The propensity of each reaction is calculated as follows:

• propensity = k * ^[r(1)]C_s(1) * ^[r(2)]C_s(2) * ... * ^[r(n)]C_s(n)

where:

• the modifier is k
  • Either a simple rate constant or a complex equation
• the reactant species for this reaction are named r(1), r(2), ..., r(n)
  • Technically, species with a location, e.g. a contained M
• the stoichiometries of each reactant are s(1), s(2), ..., s(n)
  • Positive, whole numbers
• the concentrations of each species in the relevant space are [r(1)], [r(2)], ..., [r(n)]
  • Positive, whole numbers
• ⁿC_r is the combinatoric function n! / r!(n - r)!, the number of ways of choosing r objects from a collection of n identical objects without replacement

If any of the reactants are not present in large enough numbers to satisfy the stoichiometry for the reaction (i.e. [r(i)] < s(i) for some i), then this equation is skipped and the propensity of this reaction in this location is set to zero.

This models the likelihood that all the reactants of the reaction will come close enough to react.

Each reaction takes place with respect to one membrane. All the reactant and product species must be located either around, contained, on, or under the same membrane.

For example, the following reaction is valid

                            a around M binds b on M modifier k => c contained M
                        

But the following reaction is invalid

                            # Invalid MDSL
                            a around M binds b on X modifier k => c contained Y
                        

If the initial tree defines more than one instance of the same membrane type, then each individual membrane of this type will have its own set of reactions. Each of these reactions with have their own propensities and might proceed at different rates.

For example, consider the following very simple signal transduction model

                            initial tree {simulation {20 cell}}
                            signal around cell binds receptor on cell modifier k => reporter contained cell
                        

In this model there are 20 different cell membranes. Each cell might have a different number of receptor species on it. In an environment containing a lot of signal, this reaction will have a higher propensity for those cells that have more receptors on them, so these cells will end up with more reporter inside them.

A common modelling pattern is reversible reactions.

In a reversible reaction, the reactants of the forward reaction are the products of the reverse reaction, and the products of the forward reaction are the reactants of the reverse reaction.

Reversible reactions are specified by using the <=> arrow symbol, rather than =>, and by including the modifier equation for the reverse reaction on the right hand side of the reaction.

For example, the following two code snippets define the same reactions:

                            a around cell binds r on cell modifier k_on => a_r_complex on cell
                            a_r_complex on cell modifier k_off => a around cell and r on cell
                        

                            a around cell binds r on cell modifier k_on <=> a_r_complex on cell modifier k_off
                        

Any reversible reaction could be written out as two non-reversible reactions, but using reversible reactions removes redundant MDSL code and thus reduces the chance of errors.

When specifying the list of reactants or products for a reaction, two different words can be used to combine them: binds and and.

These two words are synonyms; they are both provided so that the modeller can choose the word that best matches the biological meaning of each reaction.

For example, the following two code snippets define the same reactions:

                            a around cell binds r on cell modifier k_on => a_r_complex on cell
                        

                            a around cell and r on cell modifier k_on => a_r_complex on cell
                        

The stoichiometry of each species in a reaction is specified as a positive whole number. For example: 2 a around cell, 1 b on cell, 6 a contained cell.

If the number is not specified, then it defaults to 1. So a around cell is the same as 1 a around cell.

There are two keywords that can be used to specify a special type of stoichiometry: gene and catalyst. These are convenience keywords to say that this species should be used for calculating the propensity of the reaction (at stoichiometry 1), but this species should not be consumed by the reaction.

The two keywords gene and catalyst are synonyms. They are both provided so that the modeller can choose the word that best matches the biological meaning of each reaction.

Using one of these keywords is the same as including the species in both the reactants and products of the reaction. So the following two code snippets define the same reaction:

                            c contained cell binds x contained cell modifier k => c contained cell and y contained cell
                        

                            catalyst c contained cell binds x contained cell modifier k => y contained cell
                        

Some processes can be modelled by reactions that take time to complete. For example, consider the following very simple model of gene expression:

                            gene g contained cell binds promoter contained cell modifier k => enzyme contained cell
                        

In this model, increases in the concentration of promoter will cause enzyme to be expressed. However, the enzymes will enter the simulation as soon as the promoter binds to the gene. If timings are important in this model, then it may be desirable for the process of gene expression to take some time. i.e. the enzymes are ready to participate in other reactions a certain time after the promoter has bound the gene.

The delay keyword inserts a time delay between the reaction's reactants being removed from the simulation and its products being added to the simulation. During this time delay, other reactions can fire and the state of the simulation can change.

                            gene g contained cell binds promoter contained cell modifier k delay 5 => enzyme contained cell
                        

The delay is in units of simulation hours. Like the modifier, it can be a number, a parameter, or a complex equation.

Each reaction's modifier and delay both control the rate at which the reaction produces its products, but they control it in very different ways. The difference between the modifier and delay is as follows:

• The modifier controls how likely the reaction is to happen at each time point.
• The delay controls how much time elapses between the reactants being consumed and the products being produced.

The basic premise of the algorithm is as follows:

                            while simulation time has not elapsed
                                1. Calculate the propensities of each reaction in each membrane
                                2. Randomly choose the time at which the next reaction fires, using the reaction propensities
                                       Using Equation 21a in Gillespie's paper
                                       time to next reaction = (1.0 / (sum of propensities)) * ln(1.0 / (random number between 0 and 1))
                                3. Choose a reaction randomly, weighted by propensity
                                4. Fire this reaction
                                       Increase the simulation time by the time chosen in step 2 above
                                       Remove the reaction's reactants from the simulation
                                       Add the reaction's products to the simulation
                        

Reactions are randomly chosen and fired, with faster reactions more likely to be chosen over slower ones, with time advancing using a non-uniform timestep.

Our simulation algorithm is an optimised version of Gillespie's algorithm above.

For example, we do not recalculate the propensities of every reaction in every membrane at each step. We only recalculate the propensities of reactions that depend on species that were changed by the previous reaction, and we only recalculate propensities in membranes that could have had their species changed by the previous reaction.

The MDSL simulator has two features that are not present in the original Gillespie algorithm: membranes and delays.

Membranes allow different species to exist in different concentrations in different spaces in the same simulation. And they allow reactions to happen differently in different spaces. The original Gillespie algorithm simulates a single space. MDSL allows you to define as many spaces as required for your simulation, and have reactions move species between these spaces.

Delays allow a reaction to fire at a certain time, removing its reactants from the simulation, but waiting for a fixed time to elapse before adding its products into the simulation. The original Gillespie algorithm has reaction products added to the simulation instantaneously (delay 0 on every reaction). This is the default in MDSL, but in MDSL delays can be specified if desired.