Stochastic reaction networks in dynamic compartment populations.

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

Fig. 1.

Compartment population exhibiting chemical and compartmental dynamics. Compartment events alter the number of compartments in the population and, in general, also their content. Chemical events act only on the compartment contents, without changing the compartment number. (A) Schematic illustration of a simple example. This system is driven by an influx of compartments containing green molecules. The compartments can then randomly undergo binary fusion events. At the same time, the content of each compartment is subject to chemical modifications of two types: a bimolecular conversion of two green molecules into a yellow one, and a constant degradation of yellow molecules. (B–D) Output of one stochastic simulation of this specific model. (B) The marginals of the joint number distribution are shown at three time points. (C) Stochastic trajectory of the total number of compartments in the population. The positive and negative updates are caused by intake and fusion events, respectively. (D) The trajectories of the total amount of molecules in the population are affected by the chemical events and compartment influx, but not by compartment fusion.

1. Theoretical Results

A. Stochastic Compartment Populations

We define a compartment population as a collection of distinct entities, each being associated with its own molecular state. The state of compartment is described by a discrete-valued, -dimensional variable . A single state variable typically represents the copy number of a particular chemical species present in compartment , but may also be used to capture more coarse-grained compartment attributes such as cell types or vesicle identities. We consider the case where the population and their compartments are characterized by their molecular state only, while other physical properties such as location in space or shape are not taken into account. Thus, two compartments with the same molecular state are identical and indistinguishable in our formalism. Correspondingly, we can characterize the population in terms of the number distribution function , which counts the number of compartments that have content equal to (Fig. 1). The full state of the population is then given by the compartment number array with , which enumerates all compartment numbers within a single (and typically infinitely sized) structure of rank . The state can be understood as a multidimensional matrix, where the compartment number is found at index .

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 . Formally, we can describe those changes using stoichiometric equations of the formwhere the symbol denotes a compartment of content , and the nonnegative integers and are the stoichiometric reactant and product coefficients of transition . Furthermore, we define the arrays and , such that the population state changes by whenever transition occurs. We let denote the set of all considered transitions. Using these definitions, we can express the state of the population at any time aswith as a counting process that counts the number of occurrences of transition up to time and as the initial configuration of the system. Note that, for compactness, time dependencies are dropped in our notation in the following, but the reader should keep in mind that both and vary with time.

We next equip the counting processes with instantaneous rate functions for , which govern how likely each compartment transition happens within an infinitesimal interval of time . Throughout this work, we consider the rate functions to depend only on the current configuration of the population , consistent with Markovian dynamics. It can then be shown that the counting processes can be expressed as independent, time-transformed unit Poisson processes such thatEq. is known as the random time change representation (22). We can also write the stochastic evolution [] of the population state in differential form aswhere is the differential of the counting process , and takes value 1 whenever a transition of type occurs at time , and is zero otherwise.

B. Transition Classes

Eq. represents a continuous time Markov chain formalism for stochastic compartment populations whose dynamics are governed by an arbitrary set of transitions . While this representation is very general, it is rather impractical, because, in most relevant situations, the set 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 and into a single compartment with content , regardless of the specific value of and .

Formally, we define a transition class in two steps. First, we specify the general structure of a transition class by fixing the number of reactant compartments and the number of product compartments that are involved. For instance, in the case of compartment fusion, we would have and . Throughout this work, we will restrict ourselves to the case , but the following discussion holds true also for transitions involving more than two reactant or product compartments. We denote by and the particular choice of reactant and product compartment contents that define a distinguishable instance of class . Two transitions within a class are called distinguishable if they are associated with a different stoichiometry when expressed in the form []. In our settings, this practically means that, whenever or take value 2, we will count only pairs of or that are combinatorially distinct, because the ordering of the compartments is physically irrelevant. In order to formally enumerate the distinct transitions within a class, we introduce a bijective mapping that assigns to each distinguishable choice of and a unique index (and vice versa). This index refers to a specific stoichiometric equation of the form [], whose stoichiometric arrays and take entries and , with the symbol denoting a Kronecker delta. In practice, the mapping can be made explicit by enumerating all of the possible contents and without repeating indistinguishable instances. In the following, we will denote the image of by , which collects all of the transitions belonging to class .

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 and . This rate law can be defined in two parts. First, we introduce a probability per unit time for the reactant compartment(s) with content to participate in a transition of class , given the current state of the population, that is,with being a content-independent rate constant and being a nonnegative function which tunes the rate in terms of the reactant compartment content(s). Note that must be symmetric in its arguments when . The term is a population weight that takes into account all of the possible ways the current state could realize a transition involving reactant compartments with content . In this work, we consider to be a combinatorial weight that reflects the physical indistinguishability of compartments having equal content. Thus, we setwhich is akin to the mass action principle of standard stochastic reaction kinetics. Note that the binomials in Eq. take a value different from 1 only when .

Finally, the second part of the rate law is a conditional probability that describes how likely the reacting compartments of content result in product compartments of content . While this step might be generally of probabilistic nature, it can also be used to encode deterministic outcomes. For instance, coming back to the example of compartment fusion, the content of the product compartment is uniquely determined by fixing the contents of the two reactant compartments. Similarly to , we require to be symmetric in or whenever either of those involves two compartments. In summary, Eq. and the outcome distribution determine the propensity function of a specific transition of class , that is,which involves contents . The rate law [] provides a versatile definition, which allows us to equip a transition class with different physical properties, constraints, or selectivity. We emphasize that Eq. can be parameterized entirely in terms of the involved contents and , such as illustrated for the examples in Table 1. Additional information on how chemical reactions can be described as compartment transition classes can be found in .

Table 1.

Several examples of population transition classes and the structure of their rate laws

Description	Stoichiometry	Propensity function
Compartment intake	∅⇀h_I(n;y)[y]	h_I(n;y)=k_Iπ_I(y)
Compartment exit	[x]⇀h_E(n;x)∅	h_E(n;x)=k_Eg_E(x)n(x)
Binary coagulation	[x]+[x^′]⇀h_C(n;x,x^′,y)[y]	h_C(n;x,x^′,y)=k_Cg_C(x,x^′)n(x)(n(x^′)−δ_x,x^′)1+δ_x,x^′δ_y,x+x^′
Binary fragmentation	[x]⇀h_F(n;x,y,y^′)[y]+[y^′]	h_F(n;x,y,y^′)=k_Fg_F(x)n(x)π_F(y|x)δ_y^′,x−y
Chemical reaction	[x]⇀h_l(n;x,y)[y]	h_l(n;x,y)=k_lg_l(x)n(x)δ_y,x+Δx_l

We can now reformulate the general counting process model from Eq. using the concept of transition classes. In particular, we associate a counter with each transition within class . Moreover, we introduce the total class transition counter , which corresponds to the cumulative number of events associated with class that happened until time . The rate function of is given by the total propensity function , which follows from the superposition theorem for Poisson processes (23). The propensity represents the probability per unit time of any event in class occurring, given the current state . Based on this, we can rewrite Eq. aswhere is a finite set of transition classes. Whenever a transition in class occurs (), the state changes by a random state update with distribution . The jump process representation in Eq. 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 .

C. Stochastic Moment Dynamics

Eq. describes the stochastic evolution of the population state , under a considered set 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 , which reveals the average number of compartments in the population with content equal to . In principle, the dynamics of can be readily derived from Eq. 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 (14). As mentioned earlier, this problem is typically addressed either by neglecting these cross-correlations based on mean-field arguments or by resorting to approximate numerical schemes and Monte Carlo estimation.

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 , 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 can be defined aswith and as a vector of nonnegative integer exponents. The sum sets the order of the moment . For instance, if , then the moment corresponds to the total number of compartments present in the population, that is, . Similarly, moments of order 1 represent the total amount—or mass—of a particular species, and so forth. It is important to keep in mind that the compartment number distribution is stochastic, and therefore each population moment will be stochastic as well (Fig. 1 ). In the following, our goal is to derive an equation which captures the stochastic moment dynamics.

We begin by studying how a single transition of transition class affects an arbitrary population moment. Assume that, right before the transition, the population is in configuration , and consider an associated moment . When the transition happens, the moment instantaneously changes towith as the net change of moment due to transition of class . Correspondingly, we can write the differential change of any population moment in point process notationwhere, in the second line, is a random jump update with distribution , analogous to Eq. . Likewise, [] 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 related to the total compartment number . In particular, we will refer to any transition class that leaves the total number of compartments unchanged (i.e., ) as a chemical event, since it exclusively modifies the compartment contents. All other cases are referred to as compartment events, since for those transition classes.

D. Calculating Mean and Variance of the Population Moments

To effectively describe fluctuations in the population moments , we derived ordinary differential equations that capture the time evolution of their average and variance. We show, in , that the expectation of an arbitrary population moment satisfies the equationwhere denotes the expectation operator. Note that, since the moment change is a constant for each , the expectation could, in principle, be moved inside the second sum in []. In the latter form, however, the right-hand side (r.h.s.) of [] involves infinite sums of moments of the number distribution itself, which defeats the purpose of a low-dimensional description in terms of population moments. Instead, we show that, under certain conditions, the sum can be expressed again as a function of a finite number of population moments, so that, after applying the expectation operator, Eq. reduces to a self-contained system of differential equations. A sufficient condition for this to be the case is that 1) the function is a polynomial in and 2) the conditional distribution has moments which are polynomials in too (). In the present study, we will focus on systems which exhibit those two properties.

Analogously to Eq. 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 (). In combination with [], 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. 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 ).

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 and introducing a simplistic toy model defined by the four transition classeswhich are also illustrated in Fig. 2. The first and second transition classes in [] 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 . According to our terminology, these first two transitions classes are compartment events, since they affect the number of compartments in the population, whereas the last two transition classes in [] account for chemical modifications. The total propensities associated with [] are found to be , , , and . We start writing the stochastic differential equation (SDE) for the number of compartments in the form [],which is affected only by the occurrence of compartment events, while the chemical birth–death events do not alter . For the total population mass , we can use [] to findNote that the mass updates related to intake or exit events depend on the content of each specific transition, while, for birth or death events, they always take values or , respectively. We can express Eq. in the compact form [] too, which equalswhere we introduced the random variables with and with , which is a categorical distribution for the content of the compartment randomly exiting the system, found by evaluating . We can proceed to study the average trajectory for Eqs. and by using the result []. We obtainNot surprisingly, the evolution of is independent of that of because is just a birth–death process with constant rates and . Conversely, the dynamics of the expected total mass depends on the expected number of compartments. Corresponding equations for the variance of and can be derived analogously (). In summary, this leads to a system of six coupled ODEs which can be integrated numerically to compute the exact mean and variance of and , as shown in Fig. 2 . Note that no moment closure approximation is required for this system.

Fig. 2.

Expected moment dynamics of compartment number and total mass for the nested birth–death model and the coagulation–fragmentation case study. Solutions obtained from moment equations (ODEs) are compared with Monte Carlo averages from stochastic simulations (SSA). Error bars and shaded areas correspond to 1 SD above and below the mean. (A) Schematic illustration of the nested birth–death system. (B and C) Moment dynamics for parameters , , , and . The content of new compartments is Poisson distributed mean . The superimposed black lines show one stochastic realization of Eqs. and in B and C, respectively. A small section (highlighted in red) is enlarged in Insets to illustrate the stochastic jump dynamics. (D) Schematic illustration of the coagulation–fragmentation model. (E and F) Expected moment dynamics for parameters = 10, , , , and varying according to the ratios shown in the legend.

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, , for which we findwhere we introduced the dimensionless parameters and . The term corresponds to , the expected number of compartments at steady state, which follows from being a birth–death process (). The term in the square brackets corresponds to the steady-state average content per compartment, . The latter coincides with the first moment of the normalized expected number distribution , where we defined . We observe that setting in Eq. gives , since the intake distribution now matches the stationary distribution of the chemical birth–death process of rates and occurring in each compartment. In other words, the content of newly arriving compartments is already at steady state, thus preserving Poissonian content statistics.

To study variations across compartment contents, we calculated the variance-to-mean ratio of , where is obtained from moment equations as (). We findThe first contribution in [] 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 , Poissonian noise is recovered due to the agreement of chemical and intake processes.

Note that, in this simple linear case study, could be characterized also directly by taking the expectation on Eq. . We show, in , 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. reveals the fluctuations of the content across individual compartments, a similar analysis can be performed for the total population mass . This is demonstrated in , where we show the evolution of for three different values of and compare it against a static population without compartmental dynamics (). 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 , but we remark that multivariate scenarios can be handled analogously. We define a random coagulation class, where each pair of compartments is equally likely to fuse with rate (i.e., in Table 1). Instead, we introduce a mass-driven fragmentation class, where a compartment undergoes a fragmentation event with rate that is proportional to its content. For , we choose a uniform fragment distribution. The corresponding total class propensities read and (). 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 []. Considering these four transition classes (Fig. 2), the SDE for isbecause is equal to for any exit or coagulation event and for intake or fragmentation. The SDE for the total mass of the compartment population assumes an even simpler form,since coagulation and fragmentation events conserve mass. The random variables and are defined like in the previous case study. For space considerations, the derivation of the moment equations is left to . 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 , we plot the expected trajectories and fluctuations of and for different parameter settings and compare them to exact stochastic simulations. In all cases, we found very good agreement between both approaches. An analogous analysis of the expected second-order moment is provided in . An interesting feature emerging from Fig. 2 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 dynamicswhere we point out the dependency on the second-order moment . In a mean-field approximation, the coagulation term in Eq. would simplify to . Indeed, replacing with implies that , thereby neglecting fluctuations in the compartment number. Additionally, the linear correction in Eq. , 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 , and we equip each cell (i.e., each compartment) with a protein expression network, as shown in Fig. 3. Protein expression is described as a random telegraph process (35), where a binary gene variable can stochastically switch between an off state and on state . The active state promotes the production of a protein at rate . Furthermore, the protein is constantly degraded at rate . We define the two-dimensional compartment content variable . We account for a cell-to-cell communication mechanism where cells in the active state can promote inactive cells to switch on protein expression. Previous studies have described cell-to-cell communication by a stochastic diffusion mechanism that couples neighboring compartments on a lattice (10, 36) or, in the limit of fast diffusion, through a shared environment (37, 38). Similarly to the latter scenario, we consider that an active cell can interact with the same probability with any inactive cell in the population. Such communication mechanism can be captured by a bicompartmental transition classwhere we made explicit in the r.h.s. that, upon this transition, both cells are in the active state (details in ). The total class propensity for [] equalswhich reflects the fact that the global activation rate is proportional to the product of the number of active cells and the number of inactive cells (i.e., ) in the current configuration . Note that, in the considered model, all transition classes—including the communication events—are chemical transitions according to our definition, so that the number of compartments remains constant at its initial value . We further remark that, while the communication class [] only acts on the binary gene state , similar transitions can be defined within our framework to capture also exchange of molecules between compartments.

Fig. 3.

Moment dynamics and steady-state properties of the cell communication model and the stem cell system. (A) Schematic illustration of the reaction network and the cell–cell interaction mechanism in the cell communication model. The yellow square symbolizes the active gene state. The green square represents the expressed protein species. (B) Expected dynamics of the total protein mass for different values of , with 1 SD above and below the mean. Lines and shaded areas correspond to the result of moment equations, while dots and error bars were obtained from averaging stochastic simulations. The population comprises cells, of which only one is in the active state at time 0. Parameters are set to , , , and . (C) Variance-to-mean ratio of the steady-state total protein mass, obtained by moment equations and plotted as a function of . (Inset) The expected steady-state protein distribution in one cell, computed with stochastic simulations for the range of communication rates used in B. The dotted and dashed lines are analytical solutions, respectively, for no cell communication and infinitely fast communication (i.e., gene always active). (D) Schematic illustration of the stem cell model. Stem cells are indicated by an orange nucleus, and differentiated cells are indicated by a purple one. The green square represents a factor associated with cell cycle progression that induces stem cell division. The reference parameter values for E–H are set to , , , and . (E) (Upper) The accumulation-reset stochastic dynamics of in stem cells is shown for the initial transient of a single realization, starting with one stem cell. (Lower) The corresponding changes in total cell number and stem cell number. (F) Comparison of the expected dynamics for the total cell number (blue) and stem cell number (orange) obtained from moment equations (ODEs) and stochastic simulations (SSA), for two different initial conditions. (G) Dependency of the steady-state stem cell fraction on variations of some model parameters around their reference values, computed from moment equations. (H) Robust dynamics of the stem cell fraction, upon applying a perturbation at time where was suddenly downscaled by a factor of 5. In red, the expected stem cell fraction obtained from moment equations. In orange, one particular stochastic realization.

We are now interested in studying the protein expression dynamics in the cell population, as a function of the communication rate constant . To this end, we derived moment equations describing the averages and variances of the active cell number and the total amount of proteins (). In Fig. 3, we plot the expected total protein dynamics for different values of the communication rate . Similarly to the previous case studies, results obtained from exact SSA and the moment-based approach are in very good agreement with each other. Next, Fig. 3 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 . We note how the noise initially increases with , peaks around , and then soon starts declining toward Poissonian noise as the activation saturates. This can be appreciated also from Fig. 3, 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 , all cells in the population are constitutively maintained active, so that the protein levels are Poisson-distributed with mean . When, instead, , each cell evolves independently, and the steady-state distribution is in agreement with the canonical telegraph model, whose analytical solution is known (35). We performed an analysis equivalent to Fig. 3 for the expected number of active cells (). 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. 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 which stochastically increases over time at a constant rate . More precisely, we consider the propensity of cell division to be proportional to the current abundance of . Upon division, is reset to zero in both daughter cells, so that can be interpreted as a proxy for cell cycle progression. This is illustrated in Fig. 3, where the dynamics of across a lineage is followed over multiple rounds of cell division along a stochastic realization. This accumulation-reset mechanism provides a simple strategy to account for peaked and nonexponential interdivision time distributions (). Similar models of cell division have been considered in ref. 19, with the difference that the factor was considered to increase deterministically. Additionally, we introduce a negative feedback mechanism, which causes stem cells to differentiate at a rate that increases with their own abundance. For instance, such feedback could originate from mechanical cues due to cell crowding (44). To account for negative feedback in our model, we introduce a second-order compartment event, which mimics the interaction of a stem cell with the remaining stem cell population (). Finally, we assume that differentiated cells die or exit the system at a constant rate . Formally, the cell content of the considered model is described by , where the binary variable indicates whether a cell is either a stem cell () or a differentiated cell ().

Our goal is to study the dynamics and the variability of the total cell number and the stem cell number in the population. We remark that, in comparison to the previous case studies, the application of the moment equation method turns out to be more challenging for this model. This is because the total propensities of the division events depend on the second-order moment , which represents the total amount of in stem cells. In combination with the second-order feedback mechanism, this would lead to a large number of equations required to capture the full dynamics of all involved moments up to a certain order. Here we address this problem by combining the multivariate Gamma closure with a mean-field approximation, where correlations among certain population moments are neglected ().

In Fig. 3, 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 is affected by varying three different parameters of the model: the feedback rate , the rate , and the ratio , with held constant (Fig. 3). Interestingly, we find that the stem cell fraction is largely robust against changes in the feedback strength as well as the ratio of division rates . In regard to the former, while changing affects the number of stem cells present in the system (), the relative propensity between symmetric and asymmetric divisions remains unaffected, thereby preserving the total stem cell fraction. This is further illustrated in Fig. 3, which shows how the stem cell fraction returns to its set point upon perturbing . Instead, considering variations of , the robustness of seems to originate from the fact that the rates of symmetric divisions and feedback events compensate for each other. A more detailed analysis of the principles underlying the robustness properties of such models shall be performed in future work. This last application shows that, even though approximate, the moment equation approach provides valuable insights into the collective dynamics of cell populations.

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.