Lineage grammars: describing, simulating and analyzing population dynamics

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 environment-dependent Stochastic Tree Grammars (eSTG). It is comprised of a set of stochastic tree grammar transition rules, which are context-free 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, Predator-Prey 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 real-world 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 life-science 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 stem-cell differentiation [18]:

$$\begin{array}{c}\hfill \mathit{SC}\stackrel{0.1\frac{\mathit{\text{events}}}{\mathit{\text{day}}}}{\to }{\left\{\mathit{SC},\mathit{SC}\right\}}_{0.5}|{\left\{\mathit{Diff},\mathit{Diff}\right\}}_{0.5}\hfill \\ \hfill \mathit{Diff}\stackrel{1\frac{\mathit{\text{event}}}{\mathit{\text{day}}}}{\to }{\left\{\mathit{Diff},\mathit{Diff}\right\}}_{0.49}|{\left\{\varphi \right\}}_{0.51}\hfill \end{array}$$

In this example, SC (stem cells) divide symmetrically 0.1 times per day, while self-renewing 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:

$$\begin{array}{c}\hfill \mathit{SC}\stackrel{10\mathit{days}}{\to }{\left\{\mathit{SC},\mathit{SC}\right\}}_{0.5}|{\left\{\mathit{Diff},\mathit{Diff}\right\}}_{0.5}\hfill \\ \hfill \mathit{Diff}\stackrel{1\mathit{day}}{\to }{\left\{\mathit{Diff},\mathit{Diff}\right\}}_{0.49}|{\left\{\varphi \right\}}_{0.51}\hfill \end{array}$$

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).

Figure 1

An example of the stem cell differentiation program execution. The program was executed up to simulation time 100 days. (A) Schematic representation of the eSTG rules (without rates and probabilities). (B) Population size over time of a specific execution. (C) Cell lineage tree of a specific execution (only one cell lineage tree out of the originating SCs and Diffs is shown). (D) Average population size over time (calculated from 1000 stochastic executions). (E) Clone size distribution, which is the final population size derived from each initiating individual (calculated from 1000 stochastic executions).

Following is another example of an eSTG program for the Luria–Delbrück Model [19]:

$$\begin{array}{c}\hfill \mathit{WT}\stackrel{1}{\to }{\left\{\mathit{WT},\mathit{WT}\right\}}_{0.99}|{\left\{\mathit{WT},\mathit{MUT}\right\}}_{0.01}\hfill \\ \hfill \mathit{MUT}\stackrel{1}{\to }{\left\{\mathit{MUT},\mathit{MUT}\right\}}_1\hfill \end{array}$$

In this model, wild-type 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).

Figure 2

An example of the Luria–Delbrück program execution. The program was executed from 1 WT to 100 cells. (A) Schematic representation of the eSTG rules. (B) Typical lineage tree execution where mutations do not occur early. (C) Rare lineage tree execution where a mutation occurs early. (D) Average population size over time (calculated from 1000 stochastic executions). (E) Clone size distribution (calculated from 1000 stochastic executions). In the rare events where the mutation happens early in the lineage, the clone size of the mutated population is large.

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 stem-cell differentiation scenario:

$$\begin{array}{ll}\mathit{SC}\left(\overrightarrow{\mathit{MS}}={\overrightarrow{x}}_{\mathit{MS}}\right)\stackrel{10\mathit{days}}{\to }& \left\{\mathit{SC}\left(\overrightarrow{\mathit{MS}}=f_{\mathit{MS}}\left({\overrightarrow{x}}_{\mathit{MS}}\right)\right),\right.\\ {\left.\mathit{SC}\left(\overrightarrow{\mathit{MS}}=f_{\mathit{MS}}\left({\overrightarrow{x}}_{\mathit{MS}}\right)\right)\right\}}_{0.5}|\\ \left\{\mathit{Diff}\left(\overrightarrow{\mathit{MS}}=f_{\mathit{MS}}\left({\overrightarrow{x}}_{\mathit{MS}}\right)\right),\right.\\ {\left.\mathit{Diff}(\left(\overrightarrow{\mathit{MS}}=f_{\mathit{MS}}\left({\overrightarrow{x}}_{\mathit{MS}}\right)\right)\right\}}_{0.5}\\ \mathit{Diff}\left(\overrightarrow{\mathit{MS}}={\overrightarrow{x}}_{\mathit{MS}}\right)\stackrel{1\mathit{day}}{\to }& \left\{\mathit{Diff}\left(\overrightarrow{\mathit{MS}}=f_{\mathit{MS}}\left({\overrightarrow{x}}_{\mathit{MS}}\right)\right),\right.\\ {\left.\mathit{Diff}\left(\overrightarrow{\mathit{MS}}=f_{\mathit{MS}}\left({\overrightarrow{x}}_{\mathit{MS}}\right)\right)\right\}}_{0.49}|{\left\{\varphi \right\}}_{0.51}\end{array}$$

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]:

$$f_{\mathit{MS}}\left(x\right)=\left\{\begin{array}{l}\begin{array}{l}x+1\mathit{\text{with}}\mathit{\text{probability}}\frac{p}{2}\hfill \\ x-1\mathit{\text{with}}\mathit{\text{probability}}\frac{p}{2}\hfill \end{array}\hfill \\ x\mathit{\text{otherwise}}\hfill \end{array}\right.$$

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:

$$\begin{array}{l}\begin{array}{l}\mathit{SC}\stackrel{10\mathit{days}}{\to }{\left\{\mathit{SC},\mathit{SC}\right\}}_{0.5}|{\left\{\mathit{Diff}\left(\mathit{Gen}=1\right),\mathit{Diff}\left(\mathit{Gen}=1\right)\right\}}_{0.5}\end{array}\\ \begin{array}{ll}\mathit{Diff}\left(\mathit{Gen}=x\right)\stackrel{1\mathit{day}}{\to }& \left\{\mathit{Diff}\left(\mathit{Gen}=x+1\right),\right.\\ {\left.\mathit{Diff}\left(\mathit{Gen}=x+1\right)\right\}}_{0.49}|{\left\{\varphi \right\}}_{0.51}\end{array}\end{array}$$

Figure 3 shows various distribution statistics of the internal state Gen over the population at different time points.

Figure 3

An example of generation counter internal state. Each species of the type Diff holds an internal state called Gen which holds the number of cell divisions since the differentiation event. The histogram of the Gen values over the entire population can be calculated at different time points (e.g. after 10, 50 and 100 days, shown in (A), (B) and (C) respectively).

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:

$$A\stackrel{r}{\to }{\left\{A,A\right\}}_p|{\left\{\varphi \right\}}_{1-p}$$

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]:

$$\frac{\mathit{dN}}{\mathit{dt}}=\mathit{rN}\left(1-\frac{N}{K}\right)$$

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)

$$\frac{\mathit{dA}}{\mathit{dt}}=\mathit{Arp}-\mathit{Ar}\left(1-p\right)=\mathit{rN}\left(1-\frac{N}{K}\right)$$

For simplicity, r in the eSTG rule is the same as the r in the logistic model.

We then get:

$$\begin{array}{c}\hfill 2p-1=1-\frac{N}{K}\hfill \\ \hfill p=1-\frac{N}{2K}\hfill \end{array}$$

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):

Figure 4

An example of dynamic population growth. An example of a simple proliferation with fate probabilities and rates that are functions of the population size. (A) Schematic representation of the eSTG rules. (B) Population size over time of a logistic growth starting from a single instance. (C) The corresponding lineage representation of the specific execution. (D) Population size over time of a production-removal growth starting from a single instance. (E) The corresponding lineage representation of the specific execution.

$$\begin{array}{l}A\stackrel{r}{\to }{\left\{A,A\right\}}_p|{\left\{\varphi \right\}}_{1-p}\\ p=1-\frac{A}{200}\end{array}$$

In a different scenario, the growth is also regulated by the rate but is leading to the same steady state. Using the following production-removal equation [26]:

$$\frac{\mathit{dN}}{\mathit{dt}}=\beta -\mathit{\alpha N}$$

we can model the dynamics using the parametric eSTG by solving:

$$\frac{\mathit{dA}}{\mathit{dt}}=\mathit{Arp}-\mathit{Ar}\left(1-p\right)=\beta -\mathit{\alpha N}$$

$$r=\frac{\beta -\mathit{\alpha A}}{\left(2p-1\right)A}$$

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:

$$\begin{array}{l}A\stackrel{r}{\to }{\left\{A,A\right\}}_p|{\left\{\varphi \right\}}_{1-p}\\ r=\frac{100-A}{\left(2p-1\right)A}\\ p=\left\{\begin{array}{l}1\mathit{if}A\le 100\hfill \\ 0\mathit{if}A>100\hfill \end{array}\right.\end{array}$$

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, stem-cells are quickly amplified using self-replicating symmetric divisions, and after reaching the target size they differentiate asymmetrically into stem-cells and differentiated cells. We can describe this scenario using the following rules:

$$\begin{array}{l}\mathit{SC}\stackrel{r_1}{\to }{\left\{\mathit{SC},\mathit{SC}\right\}}_{p_1}|{\left\{\mathit{SC},\mathit{Diff}\right\}}_{\mathit{ow}}\hfill \\ \mathit{Diff}\stackrel{r_2}{\to }{\left\{\mathit{Diff},\mathit{Diff}\right\}}_{p_2}|{\left\{\varphi \right\}}_{\mathit{ow}}\hfill \\ p_1=1\mathrm{until}\left|\mathit{SC}\right|{}_{\mathit{Time}=t}=\left|\mathit{SC}\right|{}_{\mathit{Target}}\hfill \\ p_1=0\mathrm{until}\left|\mathit{Diff}\right|{}_{\mathit{Time}=t}=\left|\mathit{Diff}\right|{}_{\mathit{Target}}\hfill \end{array}$$

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:

$$\begin{array}{l}\begin{array}{l}\mathit{SC}+\frac{\mathit{dSC}}{\mathit{dt}}=\mathit{SC}+r_1{\mathit{p}}_1\mathit{SC}=\left|\mathit{SC}\right|{}_{\mathit{Target}}\end{array}\\ \begin{array}{ll}\mathit{Diff}+\frac{\mathit{dDiff}}{\mathit{dt}}& =\mathit{Diff}+r_1\left(1-p_1\right)\mathit{SC}+r_2{p}_2\mathit{Diff}-r_2\left(1-p_2\right)\mathit{Diff}\\ =\left|\mathit{Diff}\right|{}_{\mathit{Target}}\end{array}\end{array}$$

We thus extend the program with the following:

$$\begin{array}{l}{p}_1=\frac{\left|\mathit{SC}|{}_{\mathit{Target}}-|\mathit{SC}\right|}{\left|\mathit{SC}\right|\cdot r_1}\mathit{after}\left|\mathit{Diff}\right|=\left|\mathit{Diff}\right|{}_{\mathit{Target}}\hfill \\ p_2=\frac{\left|\mathit{Diff}|{}_{\mathit{Target}}-|\mathit{Diff}|-r_1|\mathit{SC}|\left(1-p_1\right)+r_2|\mathit{Diff}\right|}{2r_2\cdot \left|\mathit{Diff}\right|}\mathit{after}\left|\mathit{Diff}\right|=\left|\mathit{Diff}\right|{}_{\mathit{Target}}\hfill \end{array}$$

Figure 5 shows simulation results of a specific execution.

Figure 5

Rules for optimal development of the crypt. Simulation results of the rules for optimal development of the crypt (see main text). The rules are executed with r₁ = 1.07, r₂ = 1, |SC|_Time = 0 = 1, |Diff|_Time = 0 = 0, |SC| _Target  = 10, |Diff| _Target  = 50 (values are taken from [27]). Shown are execution results for two time windows starting with one SC. (A) Schematic representation of the eSTG rules. (B) Population size for simulation time of 10 days. The beginning of the process is shown where the switch between p₁ = 1 and p₁ = 0 occurs at around time t = 3.6. (C) The corresponding lineage representation of the specific execution. (D) Population size for simulation time of 50 days. Shown is the homeostatic phase that occurs after |Diff|_Time = t reaches |Diff| _Target at around time t = 6. (E) The corresponding lineage representation of the specific execution. It is interesting to observe the 10 clones that are maintained by the 10 SCs.

An example from a different regime is the predator/prey model of Lotka-Volterra [8]. It describes the interaction dynamics between two species using two ODEs:

$$\begin{array}{l}\frac{\mathit{\text{dPrey}}}{\mathit{dt}}=\mathit{Prey}\left(c_1-c_2\mathit{\text{Predator}}\right)\hfill \\ \frac{\mathit{\text{dPredator}}}{\mathit{dt}}=-\mathit{\text{Predator}}\left(c_3-c_2\mathit{\text{Prey}}\right)\hfill \end{array}$$

where c _i are parameters. These equations are usually translated into the following mass action kinetic reactions:

$$\begin{array}{l}\mathit{Prey}\stackrel{c_1}{\to }2\mathit{Prey},\hfill \\ \mathit{Predator}+\mathit{Prey}\stackrel{c_2}{\to }2\mathit{Predator}\hfill \\ \mathit{Predator}\stackrel{c_3}{\to }\varphi \hfill \end{array}$$

Since eSTG has only context-free 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:

Table 1 Lotka-Volterra unimolecular representation

$$\begin{array}{l}\mathit{Prey}\stackrel{r_1}{\to }{\left\{\mathit{Prey},\mathit{Prey}\right\}}_{p_1}|{\left\{\varphi \right\}}_{\mathit{ow}}\hfill \\ \mathit{Predator}\stackrel{r_2}{\to }{\left\{\mathit{Predator},\mathit{Predator}\right\}}_{p_2}|{\left\{\varphi \right\}}_{\mathit{ow}}\hfill \\ r_1=c_1+c_2\cdot \left|\mathit{Predator}\right|\hfill \\ r_2=c_2\cdot \left|\mathit{Prey}\right|+c_3\hfill \\ p_1=\frac{c_1}{r_1}\hfill \\ p_2=\frac{c_2\cdot \left|\mathit{Prey}\right|}{r_2}\hfill \end{array}$$

Figure 6 shows an example execution of the program.

Figure 6

An example execution of the Lotka-Volterra scheme. An output example of the executed program using c₁ = 2, c₂ = 0.01, c₃ = 5, |Prey|_Time = 0 = |Predator|_Time = 0 = 900. (A) Population size as a function of time. (B) A lineage tree of one of the 900 originating preys. (C) A lineage tree of one of the 900 originating predators. Both (B) and (C) exhibit the characteristic bottleneck phenomenon, where most lineages get extinct.

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 self-renewal (while keeping a constant division rate). They show that a feedback control onto the probability is a much more effective strategy for steady-state 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):

$$\begin{array}{l}\mathit{SC}\stackrel{v_0}{\to }{\left\{\mathit{SC},\mathit{SC}\right\}}_{p_0}|{\left\{\mathit{INP},\mathit{INP}\right\}}_{1-p_0}\hfill \\ \mathit{INP}\stackrel{v_1}{\to }{\left\{\mathit{INP},\mathit{INP}\right\}}_{p_1}|{\left\{\mathit{ORN},\mathit{ORN}\right\}}_{1-p_1}\hfill \\ \mathit{ORN}\stackrel{d}{\to }{\left\{0\right\}}_1\hfill \end{array}$$

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.

Figure 7

Scenarios for feedback regulation. (A) Schematic representation of the eSTG rules. Left plots - Feedback regulation onto the probability, where population size of an example execution, average population size over 1000 executions and an example of a lineage tree starting from a single SC are shown (B, C, and D respectively). Execution started with 10 SC s, ${\mathit{r}}_0=0.506,{\mathit{p}}_0=0.5,{\mathit{r}}_1=1,{\mathit{p}}_{1_{\mathit{t}=0}}=0.942,\mathit{d}=0.0138,\mathit{g}=0.0449,$ simulation time: 10. Right plots - Feedback regulation onto the rate, where population size of an example execution, average population size over 1000 executions and an example of a lineage tree starting from 10 SC s are shown (E, F, and G respectively). Execution started with 10 SC s and 200 $\mathit{INP},{\mathit{r}}_0=0.128,{\mathit{p}}_0=0.5,{\mathit{r}}_1=1,{\mathit{p}}_{1_{\mathit{t}=0}}=0.495,\mathit{d}=0.0372,\mathit{h}=0.0734$, simulation time: 20 (values are taken from [28]).

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. L-Systems [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 individual-based 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 sub-lineages 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 experimentally-obtained 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 well-known 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:

$$A\stackrel{r}{\to }{\left\{S_1\right\}}_{p_1}\left|{\left\{S_2\right\}}_{p_2}\right|\dots \left|{\left\{S_{n-1}\right\}}_{p_{n-1}}\right|{\left\{S_n\right\}}_{p_n}$$

can be converted into n separated reaction rules:

$$\begin{array}{ll}\stackrel{r\cdot p_i}{A\to S_i},\hfill & \forall i=1..n\hfill \end{array}$$

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 mass-action kinetics.

• 2. maCRN: Chemical Reaction Networks with mass-action kinetics.

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

• 4. U-gCRN: 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 mass-action law produces a system of ODEs with a special structure, here called an maODE system. In an maODE system each right-hand-side 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 non-zero 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 U-gCRN 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 U-gCRN 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·\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·\left[A\right]·\left[B\right]}{\to }C$ or to two U-gCRN reactions $A\stackrel{r\cdot \left[A\right]\cdot \left[B\right]}{\to }\mathit{C}$ and $B\stackrel{r·\left[A\right]·\left[B\right]}{\to }0$.

The family of population dynamics specifications that can be described using basic eSTG is equivalent to U-gCRN. A U-gCRN 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}{\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 U-gCRN 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 U-gCRN 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 left-hand 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₁, …, 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₀ = 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₀, from paths in ℕ* to label nodes in S₀, whose domain is non-empty and prefix-closed (that is, L(π₁, π₂) defined ⇒ L(π₁) 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 non-empty set of paths and each node has a "unique label" which is the path π that leads to it. A root-only 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}.\mathit{if}x=\mathit{nil}\mathit{then}0\mathit{else}\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)=\hfill \\ \hfill \mathit{\lambda x}.\mathit{if}x=\mathit{nil}\mathit{then}A\mathit{else}\mathit{if}x=1,\pi \mathit{then}{L}_1\left(\pi \right)\dots \mathit{else}\mathit{if}x\hfill \\ \hfill =n,\pi \mathit{then}{L}_n\left(\pi \right)\mathit{else}\mathit{undef}\hfill \end{array}$$

  where for n = 0, A = A() is a "live" leaf.

• 3. The leaf-extension operator L, π, A ⊲ (B ₁, … B _n ), which is defined if π is an A -labeled live leaf in L(L(π) = A ≠ 0), and n > 0, and B ₁ … B _n  ∈ S ₀:

  $$\begin{array}{l}L,\pi ,A⊲\left(B1,\dots \mathit{Bn}\right)=\mathit{\lambda x}.\mathit{if}x=\pi ,1\mathit{then}{B}_1,\dots ,\mathit{else}\mathit{if}x\hfill \\ =\pi ,n\mathit{then}{B}_n,\mathit{else}L\left(x\right)\hfill \end{array}$$

For example, by the above definitions a tree with root C and with n children B₁, …, B _n which are all leaves can be written as the expression C(B₁(), …, 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₁, …, B _n ) represents the tree C() where the leaf C is extended into a node with children B₁, …, B _n ; this is then the same as the tree C(B₁(), …, 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{ℝ}}^{+}$is a function from finite tuples of lineage trees to non-negative reals, with operators:

$$d\left(r,\left(L_1,\dots ,L_n\right)\right)=\mathit{\lambda x}.\mathit{if}x=\left(L_1,\dots ,L_n\right)\mathit{then}r\mathit{else}0$$

the singleton measure, which measures $\left(L_1,\dots ,L_n\right)\in {\mathbb{L}}^n$ as r and everything else as 0;

$$M_1+\dots +M_m=\mathit{\lambda x}.M_1\left(x\right)+\dots +M_m\left(x\right)$$

the sum measure, with m > 0;

$$\begin{array}{ll}L,\pi ,B⊲M& =\mathit{\lambda x}.\mathit{if}x\\ =L,\pi ,B⊲\left(L_1,\dots ,L_n\right)\mathit{then}M\left(L_1,\dots ,L_n\right)\mathit{else}0\end{array}$$

the leaf-extension measure, where π, B is a leaf in L; this is a function in ${\mathbb{L}}^1\to {\mathrm{ℝ}}^{+}$. This is the measure such that any extended tree of the form L, π, B ⊲ (L_1, …, L _n ) for some L₁, …, L _n receives the measure M(L₁, …, 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:

$$\begin{array}{c}\hfill \mathit{if}L\mathit{is}a\mathit{lineage}\mathit{tree}\mathit{and}\pi ,B\mathit{is}a\mathit{leaf}\mathit{in}L\hfill \\ \hfill \mathit{and}B\stackrel{r}{\to }{\left\{M_1\right\}}_{{\mathrm{p}}_1}|\dots |{\left\{M_m\right\}}_{{\mathrm{p}}_m}\mathit{is}a\mathit{reaction}\mathit{in}\mathit{S}\left(\mathit{the}\mathit{only}\mathit{one}\mathit{for}B\right)\hfill \\ \hfill \mathit{then}L\to L,\pi ,B⊲\left(d\left(r\cdot p_1,M_1\right)+\dots +d\left(r\cdot p_m,M_m\right)\right)\hfill \end{array}$$

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.