 Methodology article
 Open Access
 Published:
Lineage grammars: describing, simulating and analyzing population dynamics
BMC Bioinformaticsvolume 15, Article number: 249 (2014)
Abstract
Background
Precise description of the dynamics of biological processes would enable the mathematical analysis and computational simulation of complex biological phenomena. Languages such as Chemical Reaction Networks and Process Algebras cater for the detailed description of interactions among individuals and for the simulation and analysis of ensuing behaviors of populations. However, often knowledge of such interactions is lacking or not available. Yet complete oblivion to the environment would make the description of any biological process vacuous. Here we present a language for describing population dynamics that abstracts away detailed interaction among individuals, yet captures in broad terms the effect of the changing environment, based on environmentdependent Stochastic Tree Grammars (eSTG). It is comprised of a set of stochastic tree grammar transition rules, which are contextfree and as such abstract away specific interactions among individuals. Transition rule probabilities and rates, however, can depend on global parameters such as population size, generation count, and elapsed time.
Results
We show that eSTGs conveniently describe population dynamics at multiple levels including cellular dynamics, tissue development and niches of organisms. Notably, we show the utilization of eSTG for cases in which the dynamics is regulated by environmental factors, which affect the fate and rate of decisions of the different species. eSTGs are lineage grammars, in the sense that execution of an eSTG program generates the corresponding lineage trees, which can be used to analyze the evolutionary and developmental history of the biological system under investigation. These lineage trees contain a representation of the entire events history of the system, including the dynamics that led to the existing as well as to the extinct individuals.
Conclusions
We conclude that our suggested formalism can be used to easily specify, simulate and analyze complex biological systems, and supports modular description of local biological dynamics that can be later used as “black boxes” in a larger scope, thus enabling a gradual and hierarchical definition and simulation of complex biological systems. The simple, yet robust formalism enables to target a broad class of stochastic dynamic behaviors, especially those that can be modeled using global environmental feedback regulation rather than direct interaction between individuals.
Background
In recent years there has been a great interest in modeling and simulating various aspects of population dynamics in biological and ecological systems [1–4]. The increasing computational resources along with a deeper understanding of biological and ecological phenomena have led to the development of many languages for describing, analyzing and simulating concurrent stochastic processes. Many such languages specify Markovian dynamics and differ by level of abstraction, ease and complexity of the description and execution efficiency [5]. Two widely used formalisms are based on Chemical Reaction Networks (CRN) [6] and stochastic Process Algebras (PA) [7].
CRNs were originally used to describe chemical systems. A CRN description consists of a finite set of reactions acting on a finite number of species. Each reaction specifies the identity and stoichiometry of the reactants and products along with a rate constant. Many processes can be described using CRNs, for example, PredatorPrey models [8], Cellular cascade pathways [9], Cancer progression [10], Epidemics dynamics [11], and many others [1]. Each of these processes consists of a continuous interaction between individual species (the reactants) that occurs at a certain rate and produces a group of other individuals (the products, which may be empty) that can be of the same (autocatalytic) or of different type. The description of dynamical systems using CRN is relatively simple and can be used both for analytical solving and simulations. However, this approach neglects biological aspects of the described systems by treating each object (reactant or product) as a simple entity, which ignores its environmental context and structure. For example, many molecular objects maintain their overall identity while changing in specific attributes, such as chemical modification or location. When using a CRN abstraction such molecules cannot retain identity while changing state.
PAs, on the other hand, are a family of mathematical formalisms that were originally developed to model concurrent computer systems. They enable the abstraction and specification of communication and synchronization between a collection of processes by passing messages between them. One of the most well studied PA is the π calculus, which has been shown to be very useful in describing a range of biological systems [7, 12]. The language consists of processes that are mapped to realworld objects, and channels, which are mapped to communications and interactions between the different objects. A unique feature of the π calculus allows to dynamically communicate new channels between the processes (this is termed mobility), which enables the objects to keep their identity while changing their internal states or interactions with other objects. This feature is more compatible with real biological and ecological scenarios and fits well to the way we think and observe these processes. It also allows one to abstract and specify the dynamics in a more accurate fashion. It has also been shown that this abstraction can be treated as an executable computer program, allowing to stochastically simulate any specified model [13].
Many tools have been developed in order to allow and simplify the use of mathematical modeling for the lifescience community, and each one has its strengths and weaknesses [14–16]. There is no single formalism that has all the required features and choosing the appropriate one depends on the specific goals and resources of the modeler. Our goal in this work is to develop and formulate a simpler and practical tool for modeling and simulating the behavior and interaction of populations. We do so by extending the notion of Stochastic Tree Grammar (STG) [17] by incorporating both rates and probabilities to the transition rules. These can be dynamically updated by defining them as functions of the system’s state, which includes global values such as current population size, generation count or elapsed time. In addition, we extend the system by allowing each individual to hold its own internal states which can change through inheritance. We later discuss implementation of stochastic simulation and the relation to Ordinary Differential Equations (ODE).
A prominent feature of the language is that it enables to stochastically produce possible lineage trees corresponding to single executions. These lineage trees contain a representation of the entire events history of the process, including the dynamics that led to the existing as well as to the extinct individuals. As opposed to standard approaches that output only the population size dynamics, our implementation also outputs the corresponding lineage trees, which can be used to analyze the evolutionary and developmental history of the process.
Recently, Vaughan et al. [16] presented the usage of CRNs as lineage grammars and used them to simulate phylogenetic trees. Although they enable to sample possible genealogies based on the defined reaction rules, they do not allow the specification and analysis of more complex behaviors such as feedback onto the dynamic rates and general inherited properties.
Throughout the paper, we demonstrate the usability of the language by presenting a wide range of examples that can be modeled and simulated using this approach. The examples show that the language can provide simple descriptions of systems from various domains. Example parameter values were taken from the literature when available or chosen arbitrarily in order to simplify the presentation.
Results and discussion
eSTG programs
Following is an example of an eSTG program for stemcell differentiation [18]:
In this example, SC (stem cells) divide symmetrically 0.1 times per day, while selfrenewing or differentiating with the same probability (50%), and Diff (differentiated cells) can once a day either proliferate (with probability 49%) or die (with probability 51%).
Alternatively, one can define an average time to event t instead of a rate, which can be translated interchangeably into a rate using $r=\frac{1}{t}$. The above rules are then written:
An execution of an eSTG program proceeds through the stochastic application of its transition rules on its state. An example execution of the program, on an initial 10 SC and 5 Diff, can be summarized by a cell lineage tree and population size graphs shown in Figure 1B and Figure 1C. In addition to single executions, eSTG can also be used for obtaining overall population statistics, for example, to calculate the average population size over time (Figure 1D) and the distribution of clone sizes (Figure 1E).
Following is another example of an eSTG program for the Luria–Delbrück Model [19]:
In this model, wildtype bacteria (WT) are randomly mutated (in the absence of selection) to form a resistant bacteria (MUT), thus the population size of mutated bacteria varies dramatically and is dependent on the timing in which the mutation has happened. Figure 2B and Figure 2C show specific executions of typical and rare lineage trees. Averaging over many executions can yield average population size (Figure 2D) and clone size distribution (Figure 2E).
Internal states
We define internal states for each species as a vector of variables that can change, either deterministically or stochastically for each individual, with every execution of a rule. Internal states can be used to model inherited attributes, such as mutations or substance accumulation, or record historical events such as the number of generations, number of symmetrical/asymmetrical divisions, or time since historical events. We thus extend the basic rules defined above to include internal states which are functions of the predecessor’s internal states. For example, extending the previous stemcell differentiation scenario:
In this example, we define a vector of n variables $\overrightarrow{\mathit{MS}}=\left(M{S}_{1},\dots M{S}_{n}\right)$, which correspond to the number of repeats in n Microsatellite (MS) loci in the DNA [20]. In every cell division, the number of MS repeats for each locus changes according to the stochastic function f_{ MS }, which can cause either a decrease or an increase of one repeat with probability p[21]:
This simulated data can be used for example to evaluate the relationship between n, the number of MS, and the accuracy of phylogenetic reconstruction based on MS lengths of the tree (see [22, 23] for details).
Another example for the use of internal states is the following program, which counts the number of generations since each differentiation event:
Figure 3 shows various distribution statistics of the internal state Gen over the population at different time points.
Other examples of internal states can be the counting of historical events (such as how many symmetric vs. asymmetric divisions a cell went through) or measuring the time since a certain event.
Probabilities and rates as functions
Population dynamics can change based on various conditions such as population size, internal or external changes, and elapsed time. A common phenomenon in population dynamics is the reaching of a homeostasis, meaning that at a certain point, the population size reaches a steady state.
A simple example is the growth of a species until reaching a target size. Consider the following parametric rule:
Without feedback regulation on the population size, a setting of $p=\frac{1}{2}$ results in an extinction with probability 1 [24]. A simple regulation scheme is the logistic model [25]:
where N is the population size, r is the growth rate and K is the target size (also termed carrying capacity). We can use the above parametric eSTG rule to model a logistic population growth by solving:
$\frac{\mathit{dA}}{\mathit{dt}}=\frac{\mathit{dN}}{\mathit{dt}}$ (we use A as the population size of the species A)
For simplicity, r in the eSTG rule is the same as the r in the logistic model.
We then get:
Figure 4B and Figure 4C show the resulting dynamics (population size and lineage tree) starting from a single A of the following program (setting K = 100):
In a different scenario, the growth is also regulated by the rate but is leading to the same steady state. Using the following productionremoval equation [26]:
we can model the dynamics using the parametric eSTG by solving:
The steady state of this system is $\frac{\beta}{\alpha}$ and for simplicity we limit p to be either 0 or 1, and set α = 1, β = 100. We thus define the following eSTG program:
Here, the rate is inversely dependent on the population size and the population is growing until reaching the steady state that is maintained by a feedback on p, which causes either a proliferation (p = 1) or death (p = 0). Figure 4D and Figure 4E show the resulting dynamics starting from a single A.
Another interesting scenario is described in [27], where an optimal development of the intestinal crypts is analysed. In the first stage, stemcells are quickly amplified using selfreplicating symmetric divisions, and after reaching the target size they differentiate asymmetrically into stemcells and differentiated cells. We can describe this scenario using the following rules:
where X_{Time = t} is the population size of species X at time t and X_{ Target } is the target population size of X. Although not described in [27], we continue the scenario with homeostasis by solving:
We thus extend the program with the following:
Figure 5 shows simulation results of a specific execution.
An example from a different regime is the predator/prey model of LotkaVolterra [8]. It describes the interaction dynamics between two species using two ODEs:
where c_{ i } are parameters. These equations are usually translated into the following mass action kinetic reactions:
Since eSTG has only contextfree transitions, we convert the second reaction into two unimolecular reactions while preserving the 2^{nd} order rate (see Methods for a general method to convert CRNs to unimolecular reactions while preserving the same underlined ODEs). The new reactions and their rates are described in Table 1. We note that although these new reactions are not identical to the original ones, they are still in agreement with the ODEs described above. The model can be described using the following parameterized eSTG program:
Figure 6 shows an example execution of the program.
The role of different feedback strategies on the control of organ and tissue growth can be investigated through the rates and probabilities of cellular decisions. Lander et al. [28] suggest two types of feedback strategies for the Olfactory Epithelium, one on the rate of division and the other on the probability of selfrenewal (while keeping a constant division rate). They show that a feedback control onto the probability is a much more effective strategy for steadystate robustness and rapid regeneration.
The two strategies can be described using the following eSTG program (SC – Stem Cells, INP  Immediate Neuronal Precursor, ORN  Olfactory Receptor Neuron):
and the two feedback strategies are implemented by updating the INP parameters.
Strategy 1: Feedback onto the probability
${p}_{1}=\frac{{p}_{1}}{1+g\cdot \left\mathit{ORN}\right}$ where g is a constant.
Strategy 2: Feedback onto the rate:
${v}_{1}=\frac{{v}_{1}}{1+h\cdot \left\mathit{ORN}\right}$ where h is a constantFigure 7 shows possible executions generated using the two suggested strategies.
Possible extensions
Compartments
In many cases the population moves stochastically between different compartments, where each compartment corresponds to a different environment and different resources. Extending the language to include compartments allows one to define the same transition rules for species from the same type but different rates and probabilities, depending on the physical location of the individual. The system’s state is then extended to include the population size in each compartment. In addition to the regular transition rules, one also needs to define rules for the migration of each species between each two compartments.
Individual’s probabilities and rates as functions
Defining probabilities and rates for each individual separately is not recommended due to heavy computational requirements when implementing such a scenario, however, an extension of the language can support such a definition. In this case we can allow the probabilities and rates of each individual to be also dependent on its internal states. This allows each individual to have a distinct stochastic value of its probabilities and transition rates. For example, we can define a more sophisticated predator/prey model where the probability of reproduction is dependent on the individual’s age (or weight) which is represented as an internal state, or define the proliferation dynamics of a cell based on its mutations (represented as internal states).
Conclusions
Stochastic simulation is a powerful tool to execute a complicated modeling system for which a closed form analytical solution is not possible. In addition, a simulation can generate a sample of representative scenarios that can be used for further analysis or as inputs to other programs. The complexity of natural phenomena requires a formal description framework which on one hand should be rich enough to capture the complexity and dynamics of the system and on the other hand will be compact and simple so it can be widely used by a broad community and could be implemented efficiently. There are many systems that are purely generative and derive their core results by ignoring interactions (e.g. LSystems [29] and branching processes [24]). Although the assumption of independence enables certain analytical techniques, it precludes the ability to model processes and lineages that evolve through complex interactions between individuals and their environment. In order to allow both generativity and interaction, systems such as PA and CRN are more suitable. As described in [3], the trend towards individualbased stochastic models carries many advantages; they are easier to construct, more intuitive and can predict richer phenomena than population level models. In addition, it is possible to deduce population level conclusions (such as the underlined ODEs, see Methods) from the stochastic model. The presented formalism does not offer a new modeling approach in the sense that eSTG programs can be translated interchangeably into other languages (see Methods). Instead, the suggested eSTG language formalism allows a simpler description and specification of complex stochastic dynamics of individual entities. As demonstrated by the host of examples provided, these may include population level feedback from the current system’s state (either population size, internal or external factors) onto the rates and probabilities of the different species. In addition, eSTG, as a lineage grammar also enables the representation and analysis of historical events including those of extinct sublineages and transitional time points. Derivation trees produced by simulations can be examined for consistency with specific biological hypotheses [22, 30], so that eSTG models can be validated or falsified on the basis of the trees that they generate.
The language can also be used as a basis for inference and learning of the system’s governing rules, described in the eSTG formalism as the transition rules and the underlying rates and probabilities as functions of the system’s state. The question of parameter inference from biological data is an active area of study [31–34]. In our context, biological knowledge inferred from experimentallyobtained trees [22, 30, 35–43] could be used in order to infer the corresponding lineage grammars [17, 44, 45]. This will allow the use of computers and computing resources in order to gain new biological insights. This is a great challenge, especially given noise and hidden variables, and is a subject of our future work. We hope that the development of theoretical models and tools, such as the one presented here, will facilitate research in this important direction.
Methods
Stochastic simulation
eSTG programs can be naturally simulated by the wellknown Gillespie stochastic simulation algorithm [6]. Gillespie's implementation uses the rates of all possible reactions and chooses stochastically the next reaction by assuming that the time to the next reaction is exponentially distributed with rate parameters corresponding to the reaction rates.
A rule of the form:
can be converted into n separated reaction rules:
and thus existing implementations of the Gillespie algorithm can be used to determine the next reaction and the time interval. Applying these rules to build the lineage tree is described in the Operational semantics section.
The code that was used to generate the examples in this paper will be made available as an open source tool and is currently under preparation for publication.
Equivalence and conversion to other languages
In this section we compare the expressiveness of eSTG to 4 other families:

1.
maODE: Ordinary Differential Equations arising from massaction kinetics.

2.
maCRN: Chemical Reaction Networks with massaction kinetics.

3.
gCRN: Chemical Reaction Networks with general rate kinetics.

4.
UgCRN: Unimolecular Chemical Reaction Networks with general rate kinetics.
An maCRN is a chemical reaction network where each reaction has an associated rate constant, and where the instantaneous rate of a reaction is determined by the product of the rate constant with the instantaneous concentrations of the reagents. It is known that an maCRN under that massaction law produces a system of ODEs with a special structure, here called an maODE system. In an maODE system each righthandside of each differential equation for species s has the form of a polynomial over the set of species, where each monomial with a negative sign has s as a factor (raised to some nonzero power). Conversely each maODE determines a canonical maCRN that has that maODE as its kinetics. Therefore there are canonical translations back and forth between maODEs and maCRNs [46].
A gCRN is instead a Chemical Reaction Network where each reaction has an associated rate function from current or past system states to changes of concentrations. The instantaneous rate of a reaction is then given immediately by its rate function without further considerations; the class of ODEs that a gCRN may generate depends on the class of rate functions that are available.
A UgCRN is a special case of a gCRN where all the reactions are unimolecular. For sufficiently powerful rate functions it is possible to have a (nominally) unimolecular reaction depend on the concentrations of other species, so that UgCRN is in fact as expressive as gCRN. For example, an maCRN reaction $A\stackrel{r}{\to}B$ can be translated to the gCRN reaction $A\stackrel{r\xb7\left[A\right]}{\to}B$ where [A] is the instantaneous concentration of A, and an maCRN reaction $A+B\stackrel{r}{\to}C$ can be translated to a gCRN reaction $A+B\stackrel{r\xb7\left[A\right]\xb7\left[B\right]}{\to}C$ or to two UgCRN reactions $A\stackrel{r\cdot \left[A\right]\cdot \left[B\right]}{\to}\mathit{C}$ and $B\stackrel{r\xb7\left[A\right]\xb7\left[B\right]}{\to}0$.
The family of population dynamics specifications that can be described using basic eSTG is equivalent to UgCRN. A UgCRN reaction $A\stackrel{r}{\to}{B}_{1}+\dots +{B}_{n}$ can be translated into an eSTG reaction $A\stackrel{r}{\to}{\left\{{B}_{1},\dots ,{B}_{n}\right\}}_{1.0}$, and conversely an eSTG reaction $A\stackrel{r\phantom{\rule{0.25em}{0ex}}}{\to}{\left\{{B}_{1,1},\dots \right\}}_{{\mathrm{p}}_{1}}\dots {\left\{{B}_{n,1},\dots \right\}}_{{\mathrm{p}}_{n}}$ can be translated into a set of UgCRN reactions $A\stackrel{r\cdot {p}_{1}}{\to}{B}_{1,1}+\cdots ,\dots ,A\stackrel{r\cdot {p}_{n}}{\to}{B}_{n,1}+\cdots $.
The UgCRN form of eSTGs implies that one must make choices in modelling: the main species that are the focus of a model, and occur in the lefthand side of productions, will be reflected in the generated lineage trees, but auxiliary species that appear only in the rate laws will not, even when those would be considered as equal in a model based on bimolecular interactions.
Operational semantics
We will start with basic definitions for the semantics of Lineage Trees. Paths in a tree are represented as finite sequences of natural numbers π = n_{1}, …, n_{ m } ∈ ℕ* (star means finite sequence, with nil as the empty sequence, and ", " as sequence concatenation). Each number n in a path represents the n^{th} child of a node, starting from the root. Nodes in a tree are labeled by an alphabet S_{0} = S ∪ {0} consisting of species in S and a distinguished symbol 0 ∉ S (the "dead" leaf).
Definition: a tree L is a partial function in ℕ* → S_{0}, from paths in ℕ* to label nodes in S_{0}, whose domain is nonempty and prefixclosed (that is, L(π_{1}, π_{2}) defined ⇒ L(π_{1}) defined).
Definition: A leaf in a tree L is a maximal path π  one such that L(π) is defined and there is no π ' ≠ nil where L(π , π') is defined. We also say that π, B is a (B labeled) leaf in L if π is a leaf in L and L(π) = B.
Definition: A lineage tree L is a tree where each path π such that L(π) = 0 is a leaf. $\mathbb{L}\subseteq {\mathrm{\mathbb{N}}}^{*}\to {S}_{0}$ is the set of such trees.
By these definitions, a tree is a nonempty set of paths and each node has a "unique label" which is the path π that leads to it. A rootonly tree is a function from nil to some species A.
Next, we use the λ calculus notation for the definitions of lineage tree operators (if f (x) = b then we write f = λx. b). We use the element undef for partially defined functions:

1.
The lineage tree with just one dead leaf:
$$0=\mathit{\lambda x}.\phantom{\rule{0.25em}{0ex}}\mathit{if}\phantom{\rule{0.25em}{0ex}}x=\mathit{nil}\phantom{\rule{0.25em}{0ex}}\mathit{then}\phantom{\rule{0.25em}{0ex}}0\phantom{\rule{0.25em}{0ex}}\mathit{else}\phantom{\rule{0.25em}{0ex}}\mathit{undef}$$ 
2.
The lineage tree with root A ∈ S and children L _{ i }, for A ≠ 0 and n > = 0:
$$\begin{array}{c}\hfill A\left({L}_{1},\dots ,{L}_{n}\right)\phantom{\rule{0.25em}{0ex}}=\hfill \\ \hfill \mathit{\lambda x}.\phantom{\rule{0.25em}{0ex}}\mathit{if}\phantom{\rule{0.25em}{0ex}}x=\mathit{nil}\phantom{\rule{0.25em}{0ex}}\mathit{then}\phantom{\rule{0.25em}{0ex}}A\phantom{\rule{0.25em}{0ex}}\mathit{else}\phantom{\rule{0.25em}{0ex}}\mathit{if}\phantom{\rule{0.25em}{0ex}}x=1,\pi \phantom{\rule{0.25em}{0ex}}\mathit{then}\phantom{\rule{0.25em}{0ex}}{L}_{1}\left(\pi \right)\dots \phantom{\rule{0.25em}{0ex}}\mathit{else}\phantom{\rule{0.25em}{0ex}}\mathit{if}\phantom{\rule{0.25em}{0ex}}x\hfill \\ \hfill =n,\pi \phantom{\rule{0.25em}{0ex}}\mathit{then}\phantom{\rule{0.25em}{0ex}}{L}_{n}\left(\pi \right)\phantom{\rule{0.25em}{0ex}}\mathit{else}\phantom{\rule{0.25em}{0ex}}\mathit{undef}\hfill \end{array}$$where for n = 0, A = A() is a "live" leaf.

3.
The leafextension operator L, π, A ⊲ (B _{1}, … B _{ n }), which is defined if π is an A labeled live leaf in L(L(π) = A ≠ 0), and n > 0, and B _{1} … B _{ n } ∈ S _{0}:
$$\begin{array}{l}L,\pi ,A\u22b2\left(B1,\dots \mathit{Bn}\right)\phantom{\rule{0.25em}{0ex}}=\mathit{\lambda x}.\phantom{\rule{0.25em}{0ex}}\mathit{if}\phantom{\rule{0.25em}{0ex}}x=\pi ,1\phantom{\rule{0.25em}{0ex}}\mathit{then}\phantom{\rule{0.25em}{0ex}}{B}_{1},\phantom{\rule{0.25em}{0ex}}\dots ,\mathit{else}\phantom{\rule{0.25em}{0ex}}\mathit{if}\phantom{\rule{0.25em}{0ex}}x\hfill \\ \phantom{\rule{7em}{0ex}}=\pi ,n\phantom{\rule{0.25em}{0ex}}\mathit{then}\phantom{\rule{0.25em}{0ex}}{B}_{n},\phantom{\rule{0.25em}{0ex}}\mathit{else}\phantom{\rule{0.25em}{0ex}}L\left(x\right)\hfill \end{array}$$
For example, by the above definitions a tree with root C and with n children B_{1}, …, B_{ n } which are all leaves can be written as the expression C(B_{1}(), …, B_{ n }()), representing a function that given the sequence nil returns the label C, given the sequence i, nil returns the label B_{ i }, and is otherwise undefined. Similarly, the expression C(), nil, C ⊲ (B_{1}, …, B_{ n }) represents the tree C() where the leaf C is extended into a node with children B_{1}, …, B_{ n }; this is then the same as the tree C(B_{1}(), …, B_{ n }()).
A collection of eSTG reactions describes a way of generating and transforming lineage trees. We now describe how each eSTG reaction transforms a lineage tree into new lineage trees. More precisely, since eSTG reactions are stochastic/probabilistic, how each reaction produces a measure of new lineage trees, where each new lineage tree is associated with its rate of occurrence.
Definition: A measure$M\in {\mathbb{L}}^{\ast}\to {\mathrm{\mathbb{R}}}^{+}$is a function from finite tuples of lineage trees to nonnegative reals, with operators:
the singleton measure, which measures $\left({L}_{1},\dots ,{L}_{n}\right)\in {\mathbb{L}}^{n}$ as r and everything else as 0;
the sum measure, with m > 0;
the leafextension measure, where π, B is a leaf in L; this is a function in ${\mathbb{L}}^{1}\to {\mathrm{\mathbb{R}}}^{+}$. This is the measure such that any extended tree of the form L, π, B ⊲ (L_{1,} …, L_{ n }) for some L_{1}, …, L_{ n } receives the measure M(L_{1}, …, L_{ n }).
For example, C(), nil, C ⊲ (d(r, (D(), E())) + d(s, F())) = d(r, C(D(), E())) + d(s, C(F())) because C(D(), E()) has the shape C(), nil, C ⊲ (D(), E()) and so it receives measure r, and C(F()) has shape C(), nil, C ⊲ (F()) and so it receives measure s.
We are now ready to define the effect of a set of eSTG reactions S on lineage trees. This is given as a reduction relation R between lineage trees and measures. We write L → M (L reduces to M) for (L, M) ∈ R, where R is defined as the smallest relation satisfying the following rule:
This rule prescribes, for example, how to carry out a simulation of a set of eSTG reactions given an initial lineage tree: at each step apply the rule above to all applicable reactions and tree leaves, sum all the measures so obtained, and sample a new lineage tree according to the resulting measure.
References
 1.
Wilkinson DJ: Stochastic modelling for systems biology. 2011, London: Chapman & Hall/CRC Press, 2
 2.
Wilkinson DJ: Stochastic modelling for quantitative description of heterogeneous biological systems. Nat Rev Genet. 2009, 10 (2): 122133. 10.1038/nrg2509.
 3.
Black AJ, McKane AJ: Stochastic formulation of ecological models and their applications. Trends Ecol Evol. 2012, 27 (6): 337345. 10.1016/j.tree.2012.01.014.
 4.
Hucka M, Finney A, Sauro HM, Bolouri H, Doyle JC, Kitano H, Arkin AP, Bornstein BJ, Bray D, CornishBowden A, Cuellar AA, Dronov S, Gilles ED, Ginkel M, Gor V, Goryanin II, Hedley WJ, Hodgman TC, Hofmeyr JH, Hunter PJ, Juty NS, Kasberger JL, Kremling A, Kummer U, Novère NL, Loew LM, Lucio D, Mendes P, Minch E, Mjolsness ED, et al: The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models. Bioinformatics. 2003, 19 (4): 524531. 10.1093/bioinformatics/btg015.
 5.
Henzinger T, Jobstmann B, Wolf V: Formalisms for Specifying Markovian Population Models. Int J Found Comput Sci. 2011, 22 (04): 823841. 10.1142/S0129054111008441.
 6.
Gillespie DT: Stochastic Simulation of Chemical Kinetics. Annu Rev Phys Chem. 2007, 58 (1): 3555.
 7.
Regev A, Silverman W, Shapiro E: Representation and simulation of biochemical processes using the πcalculus process algebra. Pac Symp Biocomput. 2001, 6: 459470.
 8.
Fujii T, Rondelez Y: Predator–prey molecular ecosystems. ACS Nano. 2012, 7 (1): 2734.
 9.
Kholodenko BN: Negative feedback and ultrasensitivity can bring about oscillations in the mitogen‒activated protein kinase cascades. Eur J Biochem. 2000, 267 (6): 15831588. 10.1046/j.14321327.2000.01197.x.
 10.
Boman BM, Wicha MS, Fields JZ, Runquist OA: Symmetric division of cancer stem cells–a key mechanism in tumor growth that should be targeted in future therapeutic approaches. Clin Pharmacol Ther. 2007, 81 (6): 893898. 10.1038/sj.clpt.6100202.
 11.
Alonso D, McKane AJ, Pascual M: Stochastic amplification in epidemics. J R Soc Interface. 2007, 4 (14): 575582. 10.1098/rsif.2006.0192.
 12.
Phillips A, Cardelli L: Efficient, Correct Simulation of Biological Processes in the Stochastic Picalculus. Computational Methods in Systems Biology. Edited by: Calder M, Gilmore S. 2007, Berlin Heidelberg: Springer, 4695: 184199. 10.1007/9783540751403_13.
 13.
Regev A, Shapiro E: The πcalculus as an abstraction for biomolecular systems. Modelling in Molecular Biology. 2004, BerlinHeidelberg: Springer, 219266.
 14.
Ghosh S, Matsuoka Y, Asai Y, Hsin KY, Kitano H: Software for systems biology: from tools to integrated platforms. Nat Rev Genet. 2011, 12 (12): 821832.
 15.
Machado D, Costa RS, Rocha M, Ferreira EC, Tidor B, Rocha I: Modeling formalisms in Systems Biology. AMB Express. 2011, 1: 4510.1186/21910855145.
 16.
Vaughan TG, Drummond AJ: A stochastic simulator of birthdeath master equations with application to phylodynamics. Mol Biol Evol. 2013, 30 (6): 14801493. 10.1093/molbev/mst057.
 17.
Gonzalez RC, Thomason MG: Syntactic pattern recognition: An introduction. 1978
 18.
Potten CS: Stem cells. 1997, London, San Diego: Academic Press
 19.
Luria SE, Delbrück M: Mutations of Bacteria from Virus Sensitivity to Virus Resistance. Genetics. 1943, 28 (6): 491511.
 20.
Weber JL, Wong C: Mutation of human short tandem repeats. Hum Mol Genet. 1993, 2 (8): 11231128. 10.1093/hmg/2.8.1123.
 21.
Valdes AM, Slatkin M, Freimer N: Allele frequencies at microsatellite loci: the stepwise mutation model revisited. Genetics. 1993, 133 (3): 737749.
 22.
Shlush LI, ChapalIlani N, Adar R, Pery N, Maruvka Y, Spiro A, Shouval R, Rowe JM, Tzukerman M, Bercovich D, Izraeli S, Marcucci G, Bloomfield CD, Zuckerman T, Skorecki K, Shapiro E: Cell lineage analysis of acute leukemia relapse uncovers the role of replicationrate heterogeneity and miscrosatellite instability. Blood. 2012, 120 (3): 60312. 10.1182/blood201110388629.
 23.
ChapalIlani N, Maruvka YE, Spiro A, Reizel Y, Adar R, Shlush LI, Shapiro E: Comparing algorithms that reconstruct cell lineage trees utilizing information on microsatellite mutations. PLoS Comput Biol. 2013, 9 (11): e100329710.1371/journal.pcbi.1003297.
 24.
Haccou P, Jagers P, Vatutin VA: Branching Processes: Variation, Growth and Extinction of Populations. 2005, Cambridge: Cambridge University Press
 25.
Vandermeer J: How populations grow: the exponential and logistic equations. Nat Educ Knowl. 2010, 1 (8): 1
 26.
Hart Y, Antebi YE, Mayo AE, Friedman N, Alon U: Design principles of cell circuits with paradoxical components. Proc Natl Acad Sci U S A. 2012, 109 (21): 83468351. 10.1073/pnas.1117475109.
 27.
Itzkovitz S, Blat IC, Jacks T, Clevers H, van Oudenaarden A: Optimality in the development of intestinal crypts. Cell. 2012, 148 (3): 608619. 10.1016/j.cell.2011.12.025.
 28.
Lander AD, Gokoffski KK, Wan FY, Nie Q, Calof AL: Cell lineages and the logic of proliferative control. PLoS Biol. 2009, 7 (1): e1510.1371/journal.pbio.1000015.
 29.
Lindenmayer A: Mathematical models for cellular interactions in development. I. Filaments with onesided inputs. J Theor Biol. 1968, 18 (3): 280299. 10.1016/00225193(68)900799.
 30.
Reizel Y, ChapalIlani N, Adar R, Itzkovitz S, Elbaz J, Maruvka YE, Segev E, Shlush LI, Dekel N, Shapiro E: Colon stem cell and crypt dynamics exposed by cell lineage reconstruction. PLoS Genet. 2011, 7 (7): e100219210.1371/journal.pgen.1002192.
 31.
Drummond AJ, Rambaut A, Shapiro B, Pybus OG: Bayesian coalescent inference of past population dynamics from molecular sequences. Mol Biol Evol. 2005, 22 (5): 11851192. 10.1093/molbev/msi103.
 32.
Ellison AM: Bayesian inference in ecology. Ecol Lett. 2004, 7 (6): 509520. 10.1111/j.14610248.2004.00603.x.
 33.
Buckland ST, Newman KB, Fernández C, Thomas L, Harwood J: Embedding population dynamics models in inference. Statistical Science. 2007, 22: 4458. 10.1214/088342306000000673.
 34.
Reinker S, Altman R, Timmer J: Parameter estimation in stochastic biochemical reactions. IEE ProcSyst Biol. 2006, 153 (4): 168178. 10.1049/ipsyb:20050105.
 35.
Carlson CA, Kas A, Kirkwood R, Hays LE, Preston BD, Salipante SJ, Horwitz MS: Decoding cell lineage from acquired mutations using arbitrary deep sequencing. Nat Methods. 2012, 9 (1): 7880.
 36.
Reizel Y, Itzkovitz S, Adar R, Elbaz J, Jinich A, ChapalIlani N, Maruvka YE, Nevo N, Marx Z, Horovitz I, Wasserstrom A, Mayo A, Shur I, Benayahu D, Skorecki K, Segal E, Dekel N, Shapiro E: Cell lineage analysis of the mammalian female germline. PLoS Genet. 2012, 8 (2): e100247710.1371/journal.pgen.1002477.
 37.
Segev E, Shefer G, Adar R, ChapalIlani N, Itzkovitz S, Horovitz I, Reizel Y, Benayahu D, Shapiro E: Musclebound primordial stem cells give rise to myofiberassociated myogenic and nonmyogenic progenitors. PLoS One. 2011, 6 (10): e2560510.1371/journal.pone.0025605.
 38.
Siegmund KD, Marjoram P, Woo YJ, Tavaré S, Shibata D: Inferring clonal expansion and cancer stem cell dynamics from DNA methylation patterns in colorectal cancers. Proc Natl Acad Sci U S A. 2009, 106 (12): 48284833. 10.1073/pnas.0810276106.
 39.
Frumkin D, Wasserstrom A, Itzkovitz S, Stern T, Harmelin A, Eilam R, Rechavi G, Shapiro E: Cell lineage analysis of a mouse tumor. Cancer Res. 2008, 68 (14): 59245931. 10.1158/00085472.CAN076216.
 40.
Salipante SJ, Thompson JM, Horwitz MS: Phylogenetic fate mapping: theoretical and experimental studies applied to the development of mouse fibroblasts. Genetics. 2008, 178 (2): 967977. 10.1534/genetics.107.081018.
 41.
Wasserstrom A, Adar R, Shefer G, Frumkin D, Itzkovitz S, Stern T, Shur I, Zangi L, Kaplan S, Harmelin A, Reisner Y, Benayahu D, Tzahor E, Segal E, Shapiro E: Reconstruction of cell lineage trees in mice. PLoS One. 2008, 3 (4): e193910.1371/journal.pone.0001939.
 42.
Salipante SJ, Horwitz MS: Phylogenetic fate mapping. Proc Natl Acad Sci U S A. 2006, 103 (14): 54485453. 10.1073/pnas.0601265103.
 43.
Navin N, Kendall J, Troge J, Andrews P, Rodgers L, McIndoo J, Cook K, Stepansky A, Levy D, Esposito D, Muthuswamy L, Krasnitz A, McCombie WR, Hicks J, Wigler M: Tumour evolution inferred by singlecell sequencing. Nature. 2011, 472 (7341): 9094. 10.1038/nature09807.
 44.
Carrasco RC, Oncina J, CaleraRubio J: Stochastic inference of regular tree languages. Mach Learning. 2001, 44 (1–2): 185197.
 45.
Sakakibara Y: Grammatical inference in bioinformatics. IEEE Trans Pattern Anal Mach Intell. 2005, 27 (7): 10511062.
 46.
Hárs V, Tóth J: On the inverse problem of reaction kinetics. Colloquia Mathematica Societatis János Bolyai, (Szeged, Hungary, 1979) Qualitative Theory of Differential Equations. Edited by: Farkas M. 1981, 363379.
Acknowledgements
This research was supported by The European Union FP7ERCAdG Foundation and by the Kenneth and Sally Leafman Appelbaum Discovery Fund. Ehud Shapiro is the Incumbent of The Harry Weinrebe Professorial Chair of Computer Science and Biology.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
AS and ES conceived the project. AS carried out all the simulations and wrote the manuscript. LC and ES provided project oversight, formal analyses, and contributed to drafting of the manuscript and final approval. All the authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Internal State
 Lineage Tree
 Operational Semantic
 Transition Rule
 Stochastic Simulation Algorithm