# Stochastic reaction networks in dynamic compartment populations

See allHide authors and affiliations

Edited by David A. Weitz, Harvard University, Cambridge, MA, and approved August 5, 2020 (received for review February 28, 2020)

## Significance

Many biochemical processes in living systems take place in compartmentalized environments, where individual compartments can interact with each other and undergo dynamic remodeling. Studying such processes through mathematical models poses formidable challenges because the underlying dynamics involve a large number of states, which evolve stochastically with time. Here we propose a mathematical framework to study stochastic biochemical networks in compartmentalized environments. We develop a generic population model, which tracks individual compartments and their molecular composition. We then show how the time evolution of this system can be studied effectively through a small number of differential equations, which track the statistics of the population. Our approach is versatile and renders an important class of biological systems computationally accessible.

## Abstract

Compartmentalization of biochemical processes underlies all biological systems, from the organelle to the tissue scale. Theoretical models to study the interplay between noisy reaction dynamics and compartmentalization are sparse, and typically very challenging to analyze computationally. Recent studies have made progress toward addressing this problem in the context of specific biological systems, but a general and sufficiently effective approach remains lacking. In this work, we propose a mathematical framework based on counting processes that allows us to study dynamic compartment populations with arbitrary interactions and internal biochemistry. We derive an efficient description of the dynamics in terms of differential equations which capture the statistics of the population. We demonstrate the relevance of our approach by analyzing models inspired by different biological processes, including subcellular compartmentalization and tissue homeostasis.

Compartmentalization is inherent to all forms of life (1). By separating biochemical processes from their surroundings, compartments serve as spatial and functional building blocks that govern biological organization at different scales. At the subcellular level, for instance, networks of vesicles collectively regulate the delivery, sorting, and breakdown of molecular cargo (2, 3). At the tissue scale, cells themselves act as functional units, each executing its internal biochemical program while interacting with the surrounding cell population and environment. These and other forms of compartmentalization have in common that an emergent behavior or function is achieved through the collective dynamics of multiple interacting compartments, which are in complex interplay with their environment as well as the biochemical processes they carry. In many biological systems, the compartments as well as their molecular contents are present in low numbers, such that random fluctuations in their dynamics become important (4, 5). Thus, investigating the dynamical properties of compartmentalized systems, especially in the presence of random fluctuations, is an important task toward understanding living systems across different scales.

From a methodological perspective, stochasticity poses formidable challenges in the analysis of biochemical processes. Even in homogenous environments, the treatment of stochastic reactions is demanding, and a large body of literature has been devoted to addressing this subject (6⇓–8). A few previous studies, however, have attempted to study stochastic reaction dynamics within compartmentalized systems. Notable examples include the work from refs. 9 and 10, where compartments correspond to fixed spatial entities (e.g., cells in a tissue), in which reactions take place and which exchange material, for instance, through diffusion. An advantage of these models is that they can be formulated as a concatenation of multiple homogeneous reaction systems and thus could be effectively analyzed using available techniques. On the downside, however, they cannot account for compartmental dynamics, for instance, due to cell division or apoptosis.

Among the few available approaches to combine reaction and compartmental dynamics, population balance equations (PBEs) are among the most prominent (11⇓–13). In this context, PBEs describe the time evolution of the number density of a compartment population due to compartment interactions or internal compartment dynamics, which may stem from chemical reactions or material exchange. Despite their popularity and numerous applications across various fields of science (14⇓–16), PBEs are most commonly found as integro-partial differential equations in mean-field form: Rather than the actual number density, they describe the expected number density in the thermodynamic limit (17). Correspondingly, information about mesoscopic fluctuations is necessarily lost.

The relevance of noise in biological systems has led to an increased interest in stochastic population balance modeling (12, 18, 19). A few recent elegant studies, for instance, show how cell proliferation can be coupled with stochastic cell internal dynamics (20, 21). Since the adoption of a stochastic number density severely complicates the mathematical treatment, results can be often achieved only by imposing tailored approximations and simplifying assumptions, or using costly Monte Carlo simulation. In summary, while stochastic population balance has generated important insights into different biological applications, a sufficiently general yet effective framework remains lacking.

The goal of the present work is to develop a versatile and efficient approach to study stochastic biochemical processes in populations of dynamically interacting compartments. In particular, we consider both the compartments and the molecules inside them as discrete objects that can undergo arbitrary stochastic events, captured by a set of stoichiometric equations. This allows us to include interactions among distinct compartments such as compartment fusion or fission as well as chemical modifications inside each compartment (Fig. 1*A*). In contrast to population balance approaches, we describe the time evolution of a finite-size population of compartments and their molecular constituents in terms of counting processes. Based on this formalism, we then show how the population dynamics can be compactly expressed in terms of population moments. The obtained moment dynamics are themselves stochastic and therefore carry information not only about the average behavior of the population but also about its variability, as opposed to mean-field models. Using moment closure approximations, we derive a set of ordinary differential equations (ODEs), which reveal means and variances of these population moments in a very efficient manner. We demonstrate our approach using several case studies inspired by biological systems of different complexity.

## 1. Theoretical Results

### A. Stochastic Compartment Populations

We define a compartment population as a collection of N distinct entities, each being associated with its own molecular state. The state of compartment i is described by a discrete-valued, D-dimensional variable *B*). The full state of the population is then given by the compartment number array

We next allow the compartment population to exhibit temporal dynamics. On the one hand, changes in the population may occur because compartments themselves undergo modifications and interact with one another. For instance, two compartments may fuse, or a compartment may exit the system. On the other hand, a compartment’s state may change due to chemical reactions among its molecules. Regardless of their specific nature, all chemical or compartment modifications can be expressed in terms of changes in the number compartment distribution

We next equip the counting processes **3** is known as the random time change representation (22). We can also write the stochastic evolution [**2**] of the population state in differential form as

### B. Transition Classes

Eq. **4** represents a continuous time Markov chain formalism for stochastic compartment populations whose dynamics are governed by an arbitrary set of transitions J. While this representation is very general, it is rather impractical, because, in most relevant situations, the set J encompasses infinitely many transitions. For instance, if two compartments of arbitrary size can fuse with each other, then an infinite number of transitions has to be introduced to model the interaction of all possible pairs of compartment contents. To address this problem, we introduce a specification of the transitions in terms of a finite set of *transition classes*, which represent generic rules by which compartments of arbitrary content can interact. In the case of compartment fusion, for instance, we could define a single transition class which transforms two compartments with content x and

Formally, we define a transition class in two steps. First, we specify the general structure of a transition class c by fixing the number of reactant compartments **1**]. In our settings, this practically means that, whenever **1**], whose stoichiometric arrays

The second step in defining a transition class is the specification of a rate law that assigns a rate to each individual transition within the class as a function of its particular compartment contents **6** take a value different from 1 only when

Finally, the second part of the rate law is a conditional probability **5** and the outcome distribution **7**] provides a versatile definition, which allows us to equip a transition class with different physical properties, constraints, or selectivity. We emphasize that Eq. **7** can be parameterized entirely in terms of the involved contents *SI Appendix*, section S1.

We can now reformulate the general counting process model from Eq. **4** using the concept of transition classes. In particular, we associate a counter **4** as**9** is analytically convenient and, moreover, entails a natural strategy to perform stochastic simulations similar to Gillespie’s stochastic simulation algorithm (SSA) (24). Further details on how to perform stochastic simulations of the considered model are provided in *SI Appendix*, section S2.

### C. Stochastic Moment Dynamics

Eq. **8** describes the stochastic evolution of the population state n, under a considered set C of transition classes. Our goal is now to analyze the statistical properties of the population. In the context of stochastic population balance, compartment populations are typically studied based on the expected number distribution **8** by taking the expectation on both sides of the equation. Unfortunately, however, this generates an infinite-dimensional system of equations which, in general, involves higher-order cross-correlations of the number distribution such as

To circumvent these problems, we use an alternative strategy to study the dynamics of the population. In particular, we focus on *population moments*, which capture summary statistics of the full population state n, such as the total number of compartments in the population or the total number of molecules of a given type. Importantly, this will allow us to effectively access fluctuations in the dynamics of the population, while bypassing the difficulties encountered when dealing with the expected number distribution.

A moment associated with the population state n can be defined as*C* and *D*). In the following, our goal is to derive an equation which captures the stochastic moment dynamics.

We begin by studying how a single transition **9**. Likewise, [**13**] is a continuous time jump process with discrete increments.

We finally remark that a useful distinction between transition classes can be made based on the moment updates *chemical event*, since it exclusively modifies the compartment contents. All other cases are referred to as *compartment events*, since

### D. Calculating Mean and Variance of the Population Moments

To effectively describe fluctuations in the population moments M, we derived ordinary differential equations that capture the time evolution of their average and variance. We show, in *SI Appendix*, section S3, that the expectation of an arbitrary population moment satisfies the equation**14**]. In the latter form, however, the right-hand side (r.h.s.) of [**14**] involves infinite sums of moments of the number distribution **14** reduces to a self-contained system of differential equations. A sufficient condition for this to be the case is that 1) the function *SI Appendix*, section S4). In the present study, we will focus on systems which exhibit those two properties.

Analogously to Eq. **14** for the expected moment dynamics, we can derive differential equations for the expectation of a squared moment using the rules of stochastic calculus for counting processes (*SI Appendix*, section S5). In combination with [**14**], this allows us to study the expected behavior of a population moment as well as its variability across different realizations. This is an important difference from mean-field approaches, in which fluctuations in the number distribution and their corresponding moments are neglected.

We finally remark that the coupled ODE system resulting from Eq. **14** will not be closed, in general, since its r.h.s. may depend on higher-order moments. This problem can be addressed using moment closure approximations, where moments above a certain order are approximated by functions of moments up to that order (25). These approximation schemes typically rely on certain assumptions on the underlying process distribution and may give more or less accurate results depending on the details of the considered system (26). In our analyses, we found the multivariate Gamma closure as proposed in ref. 27 to give accurate results, and we will adopt this choice of closure in our case studies when needed (details in *SI Appendix*, section S6).

## 2. Case Studies

We next demonstrate our framework and the moment equation approach using several case studies inspired by biological systems at different scales. All simulations have been performed using the scientific computing language Julia (28). The code is publicly available at https://github.com/zechnerlab/StochasticCompartments.

### A. Nested Birth–Death Process

We begin by considering a population of compartments with univariate content *A*. The first and second transition classes in [**15**] are, respectively, an intake transition class, where a new compartment enters the population with a Poisson-distributed content with mean-parameter λ, and a random-exit transition class, for which any compartment can leave the population with the content-independent exit rate **15**] account for chemical modifications. The total propensities associated with [**15**] are found to be **13**],**12**] to find**17** in the compact form [**13**] too, which equals**16** and **17** by using the result [**14**]. We obtain*SI Appendix*, section S7). In summary, this leads to a system of six coupled ODEs which can be integrated numerically to compute the exact mean and variance of N and *B* and *C*. Note that no moment closure approximation is required for this system.

We next utilize the derived moment equations to investigate more closely how compartmental and chemical fluctuations affect the dynamics of the population. We first focus on the expected total mass at steady state, *SI Appendix*, section S7). The term in the square brackets corresponds to the steady-state average content per compartment, **20** gives

To study variations across compartment contents, we calculated the variance-to-mean ratio *SI Appendix*, section S7). We find**21**] corresponds to Poissonian noise stemming from the chemical fluctuations inside each compartment. The second contribution is caused by compartmental fluctuations, generally leading to super-Poissonian statistics. However, when

Note that, in this simple linear case study, **8**. We show, in *SI Appendix*, section S7, how closed-form analytical solutions can be obtained for this simple system under two special parameter configurations.

We want to point out that, while Eq. **21** reveals the fluctuations of the content across individual compartments, a similar analysis can be performed for the total population mass *SI Appendix*, section S7, where we show the evolution of *SI Appendix*, Fig. S1). This analysis shows that compartmental events are typically the dominant source of total mass fluctuations in the considered model.

This simple case study serves to illustrate how the proposed framework can be used to study fluctuations in systems that exhibit both compartment and reaction dynamics. In the following, we will consider systems with more complex interactions.

### B. Stochastic Model of a Coagulation–Fragmentation Process.

Coagulation–fragmentation processes form an important class of models to describe populations of interacting components (15), which have been used to study biological phenomena at different scales, including protein clustering (29), vesicle trafficking (30⇓–32), or clone-size dynamics during development (33). Previously, these models have been analyzed mostly using mean-field approaches or forward stochastic simulation. In this case study, we will revisit coagulation–fragmentation systems using the proposed moment equation approach. For simplicity, we will consider again a univariate compartment content *SI Appendix*, section S8). Moreover, we might consider that the population can exchange compartments with an external environment. In order to account for this, we can equip our model with an intake and an exit transition class, similarly to model [**15**]. Considering these four transition classes (Fig. 2*D*), the SDE for N is*SI Appendix*, section S8. In particular, since the considered system does not exhibit closed moment dynamics, we made use of the proposed Gamma closure as mentioned earlier. In Fig. 2 *E* and *F**,* we plot the expected trajectories and fluctuations of N and *SI Appendix*, Fig. S2. An interesting feature emerging from Fig. 2*F* is that the coagulation and fragmentation rates affect the variability of the total mass, but not its average behavior, which depends only on the intake and exit parameters. This happens because a larger coagulation rate implies that the same total mass has to be shared among fewer compartments, which causes the total population mass to exhibit larger fluctuations upon occurrence of the intake and exit events. This fact could not be captured by a mean-field treatment of a coagulation–fragmentation system, since fluctuations are necessarily lost in that case. Similar considerations hold true for the expected compartment number dynamics**24** would simplify to **24**, which originates from the exact combinatorics of the possible compartment pairings, would be omitted too. Both of these approximations can lead to significant deviations when only a few compartments are present in the system. For more details on the validity of mean-field approximations in stochastic coagulation systems, the reader may refer to ref. 34.

### C. Protein Expression Dynamics in a Cell Community

In our next case study, we apply our approach to analyze a population of compartments which are chemically coupled to each other. To this end, we consider a cell population of fixed size *A*. Protein expression is described as a random telegraph process (35), where a binary gene variable can stochastically switch between an off state *SI Appendix*, section S9). The total class propensity for [**25**] equals**25**] only acts on the binary gene state

We are now interested in studying the protein expression dynamics in the cell population, as a function of the communication rate constant *SI Appendix*, section S9). In Fig. 3*B*, we plot the expected total protein dynamics *C* illustrates the variance-to-mean ratio of the total protein amount at steady state as a function of the communication rate. This result has been obtained through moment equations, by computing *C*, *Inset*, which shows the shape of the normalized expected protein distribution in a cell, computed by SSA sampling. The distributions are in agreement with the theoretical predictions for the two limiting cases of absent and infinitely fast communication. In particular, for *B* and *C* for the expected number of active cells *SI Appendix*, Fig. S3). As can be seen from this application, moment equations provide an effective means to access the statistics of a compartment population for a wide range of parameters and identify interesting dynamical regimes with little computational effort.

### D. Stem Cell Population Dynamics

In our last application, we aim to study a system involving a more complicated interplay between internal and compartmental dynamics. Specifically, we consider a model inspired by the proliferation and differentiation dynamics of stem cell populations (39, 40), as illustrated in Fig. 3*D*. While the regulatory mechanisms controlling cell differentiation are diverse and complex, we focus here on a minimal model to exemplify how our framework can be applied to problems of such kind. We consider that stem cells can undergo division events with two probabilistic outcomes: either a symmetric division, where both daughter cells are stem cells, or an asymmetric division, where one of the daughter cells differentiates. Cell cycle length distributions are typically peaked and nonexponential (41⇓–43). To account for this, we couple the division rate of a stem cell to the dynamics of an internal molecular factor *E*, where the dynamics of *accumulation-reset* mechanism provides a simple strategy to account for peaked and nonexponential interdivision time distributions (*SI Appendix*, Fig. S4). Similar models of cell division have been considered in ref. 19, with the difference that the factor *SI Appendix*, section S10). Finally, we assume that differentiated cells die or exit the system at a constant rate

Our goal is to study the dynamics and the variability of the total cell number N and the stem cell number *SI Appendix*, section S10).

In Fig. 3*F*, we plot the expected dynamics of the total number of cells and the number of stem cells starting from two different initial conditions (1 or 100 stem cells). Even though we used additional approximations, the moment dynamics are in relatively good agreement with the results obtained from stochastic simulations. Based on the moment equations, we next investigated how the steady-state stem cell fraction *G*). Interestingly, we find that the stem cell fraction is largely robust against changes in the feedback strength *SI Appendix*, Fig. S5), the relative propensity between symmetric and asymmetric divisions remains unaffected, thereby preserving the total stem cell fraction. This is further illustrated in Fig. 3*H*, which shows how the stem cell fraction returns to its set point upon perturbing

## 3. Discussion

Compartmentalization of biochemical processes is a hallmark of living systems across different scales, from organelles to cell communities. Theoretical approaches which address the interplay between compartment and reaction dynamics are therefore of great relevance. In this work, we introduced a mathematical framework to model arbitrary compartmental and biochemical dynamics in a population of interacting compartments. Our approach relies on a fully stochastic treatment and is thus suitable to investigate the effect of mesoscopic fluctuations on compartmentalized biochemical systems. We have shown how the dynamics of a compartment population can be compactly described by ordinary differential equations, which capture means and variances of certain population moments, such as the compartment number or total molecular content. Therefore, this technique provides an analytical and computational means to efficiently access the statistical properties of the population, which would be costly to obtain using Monte Carlo simulation.

While the proposed approach is fairly general, it currently has a few limitations, which are worth addressing in the future. First, it relies on the availability of suitable moment closure approximations. In all our case studies, we found the Gamma closure to give accurate results, but different closures may be required for other types of systems. Second, we have focused on propensity functions that lead to self-contained moment dynamics, and we identified two sufficient conditions for this to be the case. While this entails a large class of systems, it will be interesting to extend our approach to more general rate laws, which are beyond these two conditions.

Our approach relies on a Markovian formalism, where the rate functions depend only on the current state of the population. However, coarse-grained events, such as cell division, can exhibit strong history dependencies. We have shown how memory effects can be included within our Markovian framework by coupling compartment events to internal processes or supplementary variables (45). For instance, in our last case study, we introduced an internal timer process which controls the rate of cell division, leading to peaked interdivision time distributions as observed experimentally (41). This strategy is similar to the work of ref. 46, where the authors have shown that a chemical system with molecular memory can be mapped onto an equivalent Markovian model.

In some of the presented case studies, we have shown how our framework can be used to track additional compartment properties in addition to their molecular content. For instance, compartments can be associated with distinct types or categories, each exhibiting different dynamical features. These categories could correspond to cell types or clusters of cells within different regions in a tissue. This could be particularly relevant for studying stochasticity in developmental systems, where cells originating from the same progenitors can commit to different fates and genetic programs in a spatiotemporal context.

## Data Availability.

Derivations and additional information are provided in *SI Appendix*. The code used for simulations is available at GitHub, https://github.com/zechnerlab/StochasticCompartments.

## Acknowledgments

We thank Quentin Vagne for his helpful comments on the implementation of stochastic simulations, and André Nadler for critical comments on the manuscript. We were supported by core funding of the Max Planck Institute of Molecular Cell Biology and Genetics.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: zechner{at}mpi-cbg.de.

Author contributions: L.D. and C.Z. designed research; L.D. performed research; L.D. and C.Z. conceptualized the theoretical approach; and L.D. and C.Z. wrote the paper.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2003734117/-/DCSupplemental.

- Copyright © 2020 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

## References

- ↵
- B. Alberts et al.

- ↵
- ↵
- ↵
- H. H. McAdams,
- A. Adam

- ↵
- W. Bialek

- ↵
- ↵
- ↵
- D. Schnoerr,
- G. Sanguinetti,
- R. Grima

- ↵
- A. J. McKane,
- T. J. Newman

- ↵
- ↵
- ↵
- D. Ramkrishna,
- M. R. Singh

- ↵
- ↵
- D. Ramkrishna

- ↵
- P. L. Krapivsky,
- S. Redner,
- E. Ben-Naim

- ↵
- M. Z. Jacobson

- ↵
- D. Ramkrishna,
- J. D. Borwanker

- ↵
- ↵
- N. Totis et al.

- ↵
- ↵
- ↵
- H. Koeppl,
- D. Densmore,
- G. Setti,
- M. di Bernardo

- D. F. Anderson,
- T. G. Kurtz

- ↵
- P. K. Andersen,
- O. Borgan,
- R. D. Gill,
- N. Keiding

- ↵
- D. T. Gillespie

- ↵
- A. Singh,
- J. P. Hespanha

- ↵
- ↵
- ↵
- J. Bezanson,
- A. Edelman,
- S. Karpinski,
- V. Shah

- ↵
- T. E. Saunders

- ↵
- ↵
- Q. Vagne,
- P. Sens

- ↵
- ↵
- ↵
- H. Tanaka,
- K. Nakazawa

- ↵
- ↵
- L. Duso,
- C. Zechner

- ↵
- H. Youk,
- W. A. Lim

- ↵
- D. T. Gonzales,
- T. Y. D. Tang,
- C. Zechner

- ↵
- T. Stiehl,
- A. Marciniak-Czochra

- ↵
- ↵
- ↵
- A. Golubev

- ↵
- C. H. L. Beentjes,
- R. Perez-Carrasco,
- R. Grima

- ↵
- ↵
- D. R. Cox

- ↵
- J. Zhang,
- T. Zhou

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Applied Mathematics

- Biological Sciences
- Biophysics and Computational Biology