A local polynomial moment approximation for compartmentalised biochemical systems

Compartmentalised biochemical reactions are a ubiquitous building block of biological systems. The interplay between chemical and compartmental dynamics can drive rich and complex dynamical behaviors that are difficult to analyse mathematically -- especially in the presence of stochasticity. We have recently proposed an effective moment equation approach to study the statistical properties of compartmentalised biochemical systems. So far, however, this approach is limited to polynomial rate laws and moreover, it relies on suitable moment closure approximations, which can be difficult to find in practice. In this work we propose a systematic method to derive closed moment dynamics for compartmentalised biochemical systems. We show that for the considered class of systems, the moment equations involve expectations over functions that factorize into two parts, one depending on the molecular content of the compartments and one depending on the compartment number distribution. Our method exploits this structure and approximates each function with suitable polynomial expansions, leading to a closed system of moment equations. We demonstrate the method using three systems inspired by cell populations and organelle networks and study its accuracy across different dynamical regimes.


Introduction
Compartmentalisation of biochemical reactions is one of the most characteristic and fundamental properties of living systems.It occurs hierarchically across length scales, ranging from organs and tissues, down to cells and subcellular structures.Compartmentalisation provides functional segregation in space and time and allows for modularity in the assembly of larger scale structure and function.
Compartmentalised reaction systems can be described as populations of enclosed volumes (e.g., organelles or cells) that evolve with time.Individual compartments contain molecules that can undergo chemical reactions, while the compartments themselves can exhibit temporal dynamics as well.As an ex-ample, two compartments may fuse or divide, or new compartments may be created or exit the system as observed in organelle networks for instance [1].The combination of chemical reactions and compartment dynamics results in challenging multiscale problems, which attracted substantial interest in the past.
Population balance equations (PBEs) provide a general mathematical framework to describe a broad class of compartmentalized reaction systems [2,3].In the context of biology, they have been applied to endosomal networks [4], heterogeneous cell populations [5,6] and clone-size dynamics [7] to name a few.PBEs are typically formulated in the thermodynamic limit and as such apply to systems where fluctuations at both the chemical and compartmental level are negligible.While stochastic systems have also been considered within the population balance framework [8,3], they are typically accesible only through computationally expensive Monte Carlo simulation.
In recent years, there has been increased interest in the theoretical and computational analysis of compartmentalized reaction systems in the presence of stochasticity [9,10,11,12,13,14].In contrast to bulk chemical systems, stochasticity arises at two levels in these systems.First, chemical reactions inside individual compartments take place at random times, leading to random variations in each compartment's state [15].Second, compartment events such as intake, division or turnover exhibit stochasticity such that the population size fluctuates with time.Most existing computational approaches to analyze compartmentalized reaction systems focus on models of cell populations, where intracellular processes are combined with population-level dynamics such as cell growth, division or phenotypic selection [16,9,17,11,18].Since the number of cells is typically large in these systems, one can consider the limit of infinitely large populations, where fluctuations in the population-level dynamics become negligible.An important consequence of this limit is that the entire population can be described by a single probability density function p(x,t), which captures how the compartment state is distributed across the population at time t.The dynamics of p(x,t) can then be derived and solved using direct numerical integration [19,20] or moment approximation techniques [11].
The situation becomes more complex when the number of compartments is small, such that the limit of large populations is no longer applicable.This can be the case for organelle networks inside cells, or small cell communities at early developmental stages.When the population size is finite, the distribution of compartment states for a particular population is a stochastic object, which at a given time t depends on the state of all compartments present in the population.Calculating the probability distribution over the full population state would be exceedingly challenging due to the combinatorial explosion of states in such systems.Moreover, such high-dimensional object would be difficult to work with and analyze in practice.To address these challenges, we have recently proposed a more effective approach, which summarizes the high-dimensional population state in terms of a small number of population moments, such as the total number of compartments, or the total amount of a certain molecule across the population [10].A key difference to the moment techniques mentioned above is that the resulting moments are stochastic, which is a consequence of the finite population size.To effectively access the statistical properties of the compartment population, we derived ordinary differential equations which capture the mean and variance of the stochastic population moments.As with conventional moment-based techniques, the resulting system of equations is not necessarily closed, which we have so far addressed using ad hoc moment-closure approximations [10].
There are two main limitations of this approach.First, identifying suitable closure functions can be difficult in practice as they are typically found in a trial-and-error fashion.Second, it is currently restricted to systems with polynomial rate laws to ensure that the resulting moment equation hierarchy is selfcontained [21,10] .In the present work, we develop a systematic moment-approximation technique which aims at addressing these two issues.The approach is based on a local polynomial approximation of the underlying rate-laws, which for the considered class of systems factorizes into a product of two functions that depend on the compartment content and compartment number distribution, respectively.The resulting scheme always leads to a closed set of moment equations and is applicable also to systems involving non-polynomial rate-laws.We demonstrate our approach using several biologically-inspired models and analyze its accuracy and runtime.Our implementation of the moment-approximation technique is publicly available through a new release of our previously developed symbolic moment generator Compartor [22].

Stochastic compartment populations
To describe the dynamics of a stochastic compartment population, we make use of our previously proposed formalism [10].We define a compartment population as a collection of N distinct entities, each one being associated with a d-dimensional content variable x ∈ X ⊆ N d .Individual dimensions of x typically correspond to the copy number of a particular chemical species, but could encode also other discrete-valued features, such as compartment categories (e.g., cell types).
The state of the population can be represented by a number distribution function n : X → N that maps any possible value of the content variable x to the number of compartments n(x) having that particular content.The number distribution can be represented by a (typically infinite) multidimensional array n = (n(x)) x∈X , where the compartment number n(x) corresponds to the element of n at index x, i.e. n x = n(x).
The compartment population can change dynamically according to a set of transition classes.Transition classes can encode changes in both the molecular content of a compartment and in the population of compartments itself: for instance two molecules may react to produce a third one within one compartment, or two compartments may fuse into a single one, merging their contents.Transition classes therefore generalise the concept of reactions to compartmentalised systems.
A transition class c is defined by a set x c = {x Consider a compartment population whose dynamics follows a set of transition classes C .The time evolution of the number distribution can then be described as where n(0) is the initial state of the population at t = 0 and R c,x c ,y c (t) are counting processes with rate function h c (n(t), x c , y c ) tracking the number of transitions of type c involving the particular sets of reactant-and product compartment x c and y c up to time t.The array ∆n c (x c , y c ) = (− ∑ x∈x c δ x,s + ∑ y∈y c δ y,s ) s∈X captures how the population state n(t) changes when this particular transition happens.The process defined in (1) can be understood as an instance of a measure-valued stochastic process [23,24].Note that for compactness, we will drop the time-dependency of n(t) and all derived quantities in the following.Throughout this work, we consider rate functions of the form with k c as a rate constant, g c as a function that depends exclusively on the content of the reactant compartments and w c,x c := w c (n(x )) as a function that depends only on the population state n.Throughout this work, we consider w c,x c to be a multilinear polynomial, which arises combinatorially by counting all indistinguishable combinations of reactant compartments for a particular population state n.We refer to functions w c,x c defined in this way as mass-action-like.The function π c denotes the outcome distribution that determines how likely a particular set of product compartments is realized from a given set of reactant compartments.In the case of cell division, for examples, π c describes the partitioning of molecular content between daughter cells given the content of the mother cell right before division.Concrete examples of transition classes and their rate functions will be provided later in our case studies.
We can observe from eq. ( 2) that the transition rate functions in the considered systems are separable: they factorise into a product of a function that depends on the content of the involved compartments and another function that depends on how many such compartments are available in the population.We can therefore write where f c (x c , y c ) only depends on the content of the reactant and product compartments and w c,x c (n) depends explicitly only on the number of compartments with the given contents (x c ), but not on the content itself.

Moment dynamics
While eq. ( 1) describes the time evolution of the full population state, this equation is difficult to deal with in practice.A more effective way to analyze compartmentalized systems is provided by moment equations, which capture certain statistical properties of the compartment population [3,10,11].We define a population moment of exponent γ ∈ N d by The order of this moment is given by |γ|, which is the sum of all elements of the multi-index γ.Note that a population moment is a functional of the stochastic number distribution n and correspondingly, is itself stochastic [10].
A special population moment is the moment of order zero M 0 , which counts the total number of compartments present in the system.In the following, we will refer to this moment as N := M 0 .Similarly, a moment of order one, i.e., M e j with e j as the j-th unit vector, counts the total amount of the j-th species across all compartments.
Starting from eq. ( 1), we can now derive a stochastic differential equation which captures the time-evolution of a population moment M γ , where ∆M γ c (x c , y c ) = ∑ s s γ ∆n c,x c ,y c (s) denotes the instantaneous change of the moment M γ when the transition with compartment sets x c and y c from class c happens.In practice, it is useful to study the expected behavior of a stochastic population moment.To this end, we can take the expectation on both sides of eq. ( 5), which yields Eq. ( 6) is an ordinary differential equation that describes the time-evolution of the expected moment ⟨M γ ⟩ and we refer to it as moment equation throughout this work.
Using the separability of h c we can rewrite eq. ( 6) as Recalling eq. ( 2), we can further split the function , where the rate function g c (x c ) only depends on the content of the reactant compartments, which is then multiplied by the probability π c of obtaining specific product compartments y c .If we now carry out the summation over y c , we obtain an expectation over ∆M γ c conditionally on x c .In particular, we obtain showing that the individual summands are again separable into two functions that depend only on the compartment content and the number distribution, respectively.In the following, we refer to f γ c and w c,x c as contentand state function, respectively.Moment equations for generic products of moments (e.g., N 2 ) exhibit the same separability, such that the following results that are based on the separability property can be applied to them as well (A).

Challenges of deriving moment equations in practice
When deriving moment equations for a particular system, we can use eq.( 9).The two functions, f γ c and w c,x c , can be directly derived from the specifications of a given transition class c and for any given moment exponent γ.However, difficulties may arise depending on their particular form.These difficulties can be best explained by going through a series of representative examples.
Linear content-and state function.The simplest scenario occurs when both the content-and state functions are at most linear polynomials.In this case the resulting moment dynamics is always closed, similarly to what happens in moment equations for bulk chemical systems.As an example, let us consider a simple degradation reaction occurring within the compartments with h(n; x) = x n x and assume that we are interested in analyzing the moments N and M 1 .The state function is given by w c,x c = n x (11) and the content functions associated with N and M 1 are and respectively.The moment equations for ⟨N⟩ and ⟨M 1 ⟩ are therefore given by d⟨N⟩ dt which can be readily solved.This is because for linear contentand state functions, the moments appearing on the right-hand side of the equations cannot be of higher order than those on the left-hand side, ensuring that the resulting moment dynamics is closed.
Nonlinear content function.We next consider a more complex scenario, where the state function is still linear, while the content function is a polynomial of higher degree.As an example, we consider the case where compartments exit the system with a rate that again linearly depends on their content: with h(n; x) = x n x .The relevant content-and state functions are now given by and the equations for ⟨N⟩ and ⟨M 1 ⟩ become d⟨N⟩ dt In this case, the equation for M 1 depends on the additional moment M 2 because the content function f 1 c (x) is a second degree polynomial.Deriving the equation for M 2 yields which in turn depends on the additional moment M 3 .More generally, it turns out that the equation for any M k depends on M k+1 , leading to an infinite-dimensional system of moment equations.This phenomenon is generally known as closure problem [21,3].
Nonlinear state function.We now consider a scenario where the state function is nonlinear.This is for instance the case when two compartments fuse, i.e., In this case, the state function is which counts the number of all possible interactions of two reactant compartments over unique combinations of x, x ′ .Let us also assume that the fusion process is independent of the reactants' content.Therefore we have h Since the propensity of this transition class is independent of the content of the two reactant compartments, the content functions are The moment equations for ⟨N⟩ and ⟨M 1 ⟩ are therefore given by d⟨N⟩ dt Notice that in the first line of eq. ( 25) the sum is taken over unique combinations of x and x ′ (referred to as {x, x ′ }), while in the remaining lines of eq. ( 25) and in eq. ( 26) the sum is rewritten to go over all combinations of x and x ′ , which requires the state function to be corrected for double-counting.Notice that eq. ( 25) for ⟨N⟩ depends on n x n x ′ , making the state function a second-degree polynomial.This causes ⟨N⟩ to depend on the moment (N) 2 = ⟨NN⟩.An equation for ⟨NN⟩ can be derived using Itô's rule for counting processes (see A): This demonstrates that in this scenario the equation for any (N) k depends on (N) k+1 , leading again to an infinite set of moment equations.Note however that this closure problem is different from the one of the previous scenario as it originates from products of moments (e.g.NNN), as opposed to higher order moments such as M k .This illustrates that the compartmentalised systems considered here display two types of closure problems, which is an important difference to bulk chemical systems or compartmentalized systems in the large N limit.
We remark that for general polynomial content-and state functions, both closure problems can occur simultaneously, leading to products of multiple higher-order moments such as NM 3 or M 2,0 M 0,3 , for instance.
Non-polynomial functions.In many practical cases, the involved rate functions may not be of polynomial form.This is the case, for instance, when considering enzyme kinetics inside compartments (e.g., Michaelis-Menten or Hill functions) or compartment fusion driven by non-polynomial coagulation kernels.In this case, the functions f γ c and/or w c,x c take more general, non-polynomial forms, causing the moment dynamics to depend on terms which cannot be expressed as moments at all.In order to study such systems using moment equations, suitable polynomial approximations are necessary.

A local polynomial (moment) approximation for compartmentalised reaction systems (LPAC)
The demonstration of the closure-related challenges outlined in section 2.3 also suggests a possible remedy.In case of polynomial content-and state functions f γ c and w c,x c , their degrees dictate the order of the moments and moment products appearing on the right-hand side of (9).Therefore, constraining the degree of f γ c and w c,x c will bound both the order of the moments and the degree of moment products.If we apply the same constraints consistently to all the moments involved, this will ultimately lead to a closed system of moment equations.
To demonstrate this idea, consider maximum degrees δ c and δ s for the content-and state function, respectively.We can first derive equations for a given set of moments using eq.( 8) and put them into their separated form.The resulting content-and state functions are then approximated by polynomials of degree δ c and δ s .Taking the expectations will lead to a linear combination of expected moments.These may involve additional expected moments, for which new equations have to be derived.The content-and state functions appearing on the right-hand side of these new equations are approximated again using polynomials of order δ c and δ s .This procedure is repeated until no new expected moments appear, resulting in a closed system of equations.
The simplest way to obtain polynomial approximations of the functions f γ c and w c,x c is to use truncated Taylor series around suitable expansion points x and n.To derive explicit expressions for the approximate moment equations, we first define the following shorthand notation where x c is the set of reactant compartments of size r c corresponding to the c-th transition class and (x c ) k,i denotes the i-th chemical species of the k-th compartment in x c .The exponents α k,i and β k denote the differentiation orders associated with the content-and state functions, respectively.
The Taylor expansions of the content-and state functions can then be written as and respectively.
If we substitute the expansions, eqs.( 33) and ( 35), into the separated moment equation, eq. ( 9), we finally obtain the full explicit form of the Taylor-based moment approximation We now choose the first expasion point as x = ( xi where e i is the i-th unit vector.This can be understood as an estimate of the mean compartment content.For the second expansion point we choose the expected number distribution This is particularly important as it allows us to conveniently translate the ∑ x x γ η terms arising in the Taylor expansion (36) into functions of moments as We can now plug the chosen expansion points into (36) to obtain the final compact form:  where we used the standard multi-index notation for all operators involving tuples to arrive at the final compact form.See C for the precise definition of all multi-index operations used in our work.Note that eq. ( 40) relies on the assumption that the state function w c,x c is multilinear, as indicated by the summation domain β ≤ e, i.e. β i ≤ 1 for all i.If we now truncate eq. ( 40) such that |α| ≤ δ c and |β | ≤ δ s we obtain the final approximated moment equation for ⟨M γ ⟩.

Case studies
To test the proposed moment expansion method, we applied it to three different case studies.These examples are chosen to resemble characteristic features of biological systems such as non-linear propensities, interactions among compartments and feedback loops.We analyze the accuracy of our momentapproximation approach, by comparing it to exact Monte Carlo simulation.Moment equations were derived automatically by the Compartor [22] toolbox, which was extended to support the presented LPAC moment-expansion method.The resulting equations were solved using the DifferentialEquations.jl[25] package of the Julia programming language [26].

Binary birth-death-fusion process
We first tested our method on a simple system that features two nonlinear propensity functions.Compartments contain a single chemical species X which is synthesized with a constant rate k b .For molecular turnover, we consider an annihilation reaction 2X → / 0 with rate constant k d .Compartments enter the system at a rate k I and their initial content is distributed according to an intake distribution π I , which we consider to be a Poisson distribution with mean λ .Moreover, two compartments can fuse with each other with rate constant k F .In total, the system is given by the transition classes / 0 h I (n;y)

− −−− → [y]
(Intake) with propensity functions where x, x ′ and y are variables representing the compartments' chemical content.Here we observe that both the fusion-and annihilation transitions involve propensities that are nonlinear in state and content, respectively.Our goal is to use our moment-expansion method to study the mean and variance of the population moments N and M 1 .To this end, we require equations for the averages ⟨N⟩ and M 1 as well as the second-order moments N 2 and (M 1 ) 2 Chemical content from which variances can be obtained subsequently.However, deriving the equations for these moments leads to a closure problem due to the nonlinear propensities.We therefore employ a second-order polynomial approximation based on the Taylor series, which leads to a closed system of equations for the moments ⟨N⟩, N 2 , M 1 , (M 1 ) 2 , M 2 , (M 2 ) 2 , NM 1 , NM 2 and M 1 M 2 .For further details refer to F. Figure 1 shows that the approximate statistics of N and M 1 (dashed and dotted lines) are in close agreement with Monte Carlo simulations (solid lines and shaded areas).Quantitative error measures are reported in table 1.On our installation, generating n = 100 Monte Carlo samples took 4.7s while solving the moment equations took 6.6 • 10 −4 s, corresponding to a speedup of about 7121 times.

Shared Antithetic Integral Controller
As a second case study we apply our approximation approach to a population-level feedback control motif that we analyzed previously [12].This system is based on the antithetic integral control motif -a general biochemical control strategy that ensures robust perfect adaptation [27].In contrast to the original scheme, the considered control system -termed shared antithetic integral controller (sAIC) -acts at the multicellular level and thereby provides a means to robustly regulate population-level features such as the total amount of a produced species or the number of cells in the population.In our original work [12] we applied the sAIC to a simple cell population model and we adopt a similar model for the present case study.
Here we consider a population of cells undergoing division and apoptosis and within these cells a compartmentalised species Q is produced and degraded.This system is extended by the sAIC circuit by introducing two bulk controller species Z 1 , Z 2 and four additional transition classes.In total, the considered system is given by the following transition network with propensity functions As we have shown previously, this circuit regulates the average total number of molecules ⟨Q⟩, i.e.M 1 , to tunable set points Q * = k re f /k meas , irrespective of all other system parameters.
For further details, the reader may refer to our original work [12] and G.
Figure 2 shows a controller with setpoints Q * that change at time t = 200 and t = 400, and the total mass of the controlled species Q follows accordingly.Moreover, the statistics of the cell number and the total number of molecules are again approximated accurately by the moment-expansion method.Quantitative error estimates can be found in table 1.

Mutually repressing gene circuit in a cell population
In our last case study, we want to test our approach for a system that contains multiple chemical species that exhibit complex correlations.To this end, we take inspiration from developmental biology and consider a population of dividing cells that express two mutually repressing genes.Cell turnover is considered to be contact-dependent, where cells undergo apoptosis upon interacting with another cell from the population.The overall system is given by the transition classes [x] [x] [x] with propensity functions where x, x ′ and y are variables representing the compartments' chemical content, k D , k A , k p and k d are the rate constants for the respective transition classes, k R 1 and k R 2 are Michaelis-Menten constants and π D is the outcome distribution for cell division, where content is uniformly distributed among daughter cells.Since the two internal species mutually repress each other (eq.( 55) and (57)), we generally expect negative correlations between them.Here we want to study how the cell population dynamics impacts this correlation and to what extent our moment-approximation method is able to capture this.As a measure for correlation, we use the Pearson correlation coefficient, which can be estimated from the approximate population moments as (for details see E). Figure 3, panels (a-d), shows the dynamics of the system in a regime where cell turnover and gene expression take place with a timescale ratio of k b /k D = 200, which is roughly the magnitude at which many physiological processes take place [28].The respective quantitative error estimates can be found in table 1. Panels (e-i) show how the distribution of compartment content, together with their correlations, change over a spectrum of different timescale ratios.We speculate that the approximation (dashed lines) of moments M 1,0 , M 1,1 and ρ (panels b-d) do not match the Monte Carlo estimation (solid lines) pefectly, likely because the distributions are more disperse than in the other case studies we presented in this work (see Appendix figs. 4 and 5).Thus, while the approximation is still able to capture the qualitative dynamics of the system, this case study also reveals the potential limitations of the approach when dealing with strongly non-linear systems with disperse populations.

Discussion
Compartmentalization of biochemical reactions into compartments is a hallmark of biological systems.Studying the dynamic interplay between chemical reactions and compartmental dynamics is mathematically challenging, especially in the presence of noise.In this work we have developed LPAC, a systematic moment approximation to analyze the statistical properties of compartmentalized reaction systems.Our approach is based on our earlier work [10], but replaces ad hoc moment-closure approximations by polynomial approximations.This renders the approach applicable to a broader class of systems, including systems involving non-polynomial rate-laws.Our approach exploits the particular structure of the underlying moment equations, where individual terms can be written as expectations of products of two functions, one depending exclusively on the compartment content and one depending exclusively on the number distribution.We have further shown that these two functions give rise to two different types of closure problems, which is a key difference to bulk chemical systems as well as compartmentalized systems in the large N limit.The proposed moment-approximation technique addresses both of these closure problems and therefore always leads a closed system of moment equations.We have included this approach in a new release of our previously developed moment generator Compartor [22], which can be used to derive moment equations in a fully automated manner for a broader class of systems.
To analyze the performance of our method, we have applied it to several models exhibiting nonlinear dynamics.We generally found that the approximate moments agreed very well with those obtained from exact Monte Carlo simulations.This is especially the case for the first two case studies, in which the compartment content was narrowly distributed (Figures 1, 2, 4  and 5).However, quantitative deviations were observed in situations where the distribution of molecular content is strongly disperse (Figures 3 and 6).This is in line with the fact that a simple truncated Taylor expansion was used as a polynomial approximation, which is guaranteed to be accurate only locally around the chosen expansion point.We remark, however, that our approach could in principle be extended to other polynomial approximations.Polynomial interpolation, for instance, could be used to achieve better approximation accuracy across larger domains, if a suitable set of interpolation nodes is chosen.This in turn could achieve accurate results even in the presence of strongly disperse or multimodal distributions.Also Bernstein polynomial may provide a promising approximation technique.As shown by Lunz et al. [11], they can be numerically better behaved than interpolants, although it is currently unclear how they perform on the particular class of systems that we considered in this work.Exploring these techniques more deeply in the context of these systems is an important goal for the future as it may expand the scope of our approach to an even larger class of systems.

Supplementary Material A Equations for products of moments
As mentioned in the main text, it is possible to derive moment equations for generic products of moments.Let us consider an m-ary product of moments Using Itô's rule for counting processes [29], it is straightforward to show that the expectation of this moment satisfies the differential equation We can expand the products of all binomials in eq.(A.3) and get a sum of products of moments and ∆M where if we define By substituting the definition of a population moment into eq.(A.4) we get that each term in eq.(A.3) can be written as follows: Here x K is a tuple of compartment content symbols involved in the cross product identified by the set K , in the same way as x c represents the tuple of reactant compartments of a given transition c.
Now we can rewrite the whole moment equation for a generic product of moments as: This leads us directly to the following: Computing the contribution of a r c -compartmental transition to the equation of a product of m moments is equivalent to a "virtual" transition involving r c +m−1 compartments, with a propensity that is still separable into a content-dependent function f and a state-dependent function w.
Remark 2. This means that all the results on approximating moment equations of the form can also be applied to any arbitrary product of moments and therefore to derive the equation for the expectation of any arbitrary moment polynomial.

B Approximated moment equations
Let us now consider the separated form of a generic moment equation (9) and assume that all content and state functions have been approximated by suitable polynomials.Expanding the products leads to monomials of the form where C contains all the constant terms of the expression.Since the expectation and summation ∑ x c • • • are linear, we can bring both of them in and around each monomial.With a rearrangement of the operands and then splitting the summation we get: (B.9) Notation.Here we have used an extended notation for generalized population moments: With the term order we still refer to |γ|, while we refer to |σ | as n-order.
Remark 3. If the compartment events follow mass-action-like kinetics, then their state function w is multilinear, i.e. it is linear with respect to each single reactant compartment symbol.This means that it can only generate moments that have a maximum n-order of 1 each.
Lemma 2. If the polynomial approximating the content function f has degree δ , then for any moment M γ,σ that appears in the approximated moment equation it holds |γ| ≤ δ .
Proof.Since the degree of the polynomial is δ , then by definition the sum of the α exponents of all monomials in eq.(B.7) is bounded by δ .This means that the order of all moments generated by these monomials is also necessarily bounded by δ .
Remark 4. From the proof we also get a stronger result: the total order of the moments in any moment product is also bounded by the approximation degree δ .
Lemma 3. If the polynomial approximating the content function w has degree ε, then the total n-order of any moment product Proof.Since the degree of the polynomial approximating w is ε, then all monomials in eq.(B.7) have a total order of the β exponents that is bounded by ε.Following the derivation of eq.(B.8) we then obtain that any product of moments obtained by the approximation and truncation procedure has a total n-order bounded by ε.
Remark 5.This also means that the number of moments involved in any product appearing on the right-hand side is also bounded.In case the compartment events follow a mass-actionlike kinetics, then the number of moments must be ≤ ε.
The above results lead to the following: Theorem 1.Given an equation for any arbitrary moment or product of moments and two approximation degrees δ , ε ∈ N, the approximation of order (δ , ε) of such equation can be solved through a system of moment equations with the same approximation order and the system is closed.
Proof.Using lemmas 2 and 3 we obtain that the equation for any desired moment, approximated at order (δ , ε) only contains products of moments of total order δ and total n-order ε.Furthermore any product of moments additionally required in the system of equations can also be described with an approximated equation only containing products of moments of total order δ and total n-order ε.Since order and n-order are positive integers and bounded by δ and ε respectively, then the combinations of moments that can appear in a product of moments are finite and this proves that the system of equations is closed.
We therefore have a simple algorithm allowing to always derive a closed set of moment equations: first derive the equations for all the moments of interest, substituting all contentand state functions with suitable polynomials of limited degree; then collect all moments and moment products these equations depend on and iteratively apply the same procedure for all those we still have no equation for.This procedure will terminate and produce a closed set of moment equations.

C Multi-index notation and operators
Throughout this work, we rely on standard multi-index notation to make the equations more compact and easier to read.Multiindex notation allows common operators that accept integer indices as arguments to be extended to ordered tuples of indices.
Here we provide the definitions of multi-index operators used in our work.
A d-dimensional multi-index is defined as a d-tuple: If α, β ∈ N d 0 and x ∈ R d we can define: Partial derivatives

D Including bulk chemical species and reactions
It is possible to have systems where a compartment population, with its internal chemistry, interacts with chemistry occuring in the bulk of the outer medium.We gave an example coming from bio-engineering in section 3.2, but many other biological systems have this characteristic too.
Here we want to show how the same mathematical framework that we used can be extended in order to account for bulk chemical species and how their chemistry can be coupled to the compartmentalised system and its transitions.
Let us recall the generic 3-term form of a moment equation, from eq. ( 8): We introduce a new stochastic state vector B, which collects the copy numbers of all bulk chemical species.We can now account for the bulk chemical state in the moment equations by assuming the transition propensities also have a dependence on the bulk state: Where the bulk-dependent rate function b c would have a constant value of 1 in each transition class that does not depend on the bulk chemical state.Analogously to the presented moment equations, we can also derive differential equations for the evolution of the expected bulk state: Chemical reactions that occur in the bulk, without involving compartments and their content, are represented by transition classes that do not involve any reactant compartment, i.e. with x c = {}.In this case the corresponding terms in eq.(D.13) contain a ∆B c that is independent of x c such that the inner sum simplifies to ∆B c ⟨b c (B)⟩.In the special case where there are no chemical interactions between any bulk and compartmentalised chemical species, eq.(D.13) would further simplify to a conventional moment equation for bulk chemical systems d⟨B⟩ dt The function of the bulk state b c (B) can then be polynomially expanded around the mean state ⟨B⟩ in the same way as the state function w c,x c (n), obtaining a closed set of moment equations that account for bulk quantities and moment-bulk cross products.

E Pearson correlation coefficient in compartment populations
We use a mean-field approximation of the Pearson correlation coefficient in order to be able to compute it from the population moments that are already available in the Mutual Repression case study from section 3.3.We use the following base approximations: Plugging them into the Pearson correlation's definition we get the formula that we use in the analysis of the case study: . (E.17)

F Case study: Binary birth-death-fusion process
As detailed in section 3.1 of the main text, the system is defined by the transition classes / 0 with propensity functions The complete set of moment equations for this case study, together with the code to generate them automatically with Compartor, can be found in the corresponding Jupyter notebook in the Compartor repository on GitHub.

G Case study: Shared Antithetic Integral Controller
The sAIC case study features a bulk chemistry that is coupled to the compartmentalised system.As detailed in Supplementary Material appendix D, bulk chemistry can be integrated into the mathematical framework with little symbolic overhead.Compartor, our tool for generating the moment equations, does however not yet support this extension.We have therefore resorted to "virtualising" the sAIC bulk chemistry into a set of compartmentalised transitions that have identical dynamics.
The internal compartment state is given by the 3-species tuple (Z 1 , Z 2 , Q), where Z 1 and Z 2 are the two bulk controller chemical species and Q is the actual internal measured species.Species Z 1 and Z 2 are therefore virtually spread among the compartments in the system and their total amounts are given by the moments M 1,0,0 and M 0,1,0 respectively.The transition classes are modified in order to conserve Z 1 and Z 2 in case of compartmental events.Furthermore, the Comparison reaction is split into two transitions, in order to account for the reactant molecules to virtually belong either to the same compartment or to two different ones.
The full augmented transition network is represented by the following diagram: while the propensity functions are as follows: The complete set of moment equations for this case study, together with the code to generate them automatically with Compartor, can be found in the corresponding Jupyter notebook in the Compartor repository on GitHub.

H Case study: Mutually repressing gene circuit in a cell population
This case study features two internal chemical species, representing levels of gene expression, each inhibiting the production of the other through a repressor Michaelis-Menten type of kinetics.This internal dynamics is then imposed on a population of dividing cells, which size is kept in check by a binary cell death, i.e. a bi-compartmental transition in which one of the reactant compartments leaves the population.
The full transition network is represented by the following diagram: Figure 5: Histograms for the distributions of (a,b) the bulk controller species Z 1 and Z 2 respectively, (c) the number of cells in the system N, (d) the cellular content variable q and (e) the total mass of q in the system Q = M 1 in the Shared Antithetic Integral Controller case study.The solid black lines in each panel shows the estimated values of ⟨Z 1 ⟩, ⟨Z 2 ⟩, ⟨N⟩, ⟨Q⟩ / ⟨N⟩ and ⟨Q⟩, respectively, as computed by the approximation.
While the propensity functions are as follows: The complete set of moment equations for this case study, together with the code to generate them automatically with Compartor, can be found in the corresponding Jupyter notebook in the Compartor repository on GitHub.
reactant compartments, a set y c = {y (1) c , . . ., y (p c ) c } of product compartments and by a rate function h c (n(t), x c , y c ) which determines how likely a transition with a particular combination of reactant and product compartments occurs per unit time.Throughout this work, we consider the compartment population to exhibit Markovian dynamics such that the rate functions depend only on the current state of the population n(t).

Figure 1 :
Figure 1: Comparison of the results obtained by averaging 100 trajectories of the stochastic simulation algorithm (SSA), shown by the solid line and shaded area, with the predictions of our moment-expansion method, shown by the dashed line and dotted boundaries.The solid and dashed lines denote average values, while the shaded and dot-bordered areas are the regions within one standard deviation.Panel (a) shows the statistics of the number of compartments (N) and panel (b) shows the statistics of the total number of molecules in the system (M 1 ).

Figure 2 :
Figure 2: Comparison of the results obtained by averaging 8 SSA trajectories, shown by the solid line and shaded area, with the predictions of our momentexpansion method, shown by the dashed line and dotted boundaries.The solid and dashed lines denote average values, while the shaded and dot-bordered areas are the regions within one standard deviation.Panel (a) shows the statistics of the number of compartments (N) and panel (b) shows the statistics of the total number of molecules of different species in the system (Z 1 , Z 2 and Q), together with the controller's setpoint Q * indicated by the grey dotted line.
where x, x ′ and y are variables representing the compartments' chemical content, k D , k A , k p , k d , k re f , k act , k meas and k comp are the rate constants for the respective transition classes and π D is the outcome distribution for cell division, where content is uniformly distributed among daughter cells.The Division and Apoptosis transition classes capture the dynamics of the cell population, where individual cells can divide and turn over.Production and Degradation describe the dynamics of the cell-internal species Q. Feedback-control is achieved via the four additional transition classes: Reference and Actuation are responsible for setting the target level of Q and promoting its production, while Measurement and Comparison provide a way to sense the current value of Q and downregulate its production.

Figure 3 :
Figure 3: Comparison of the results obtained by averaging 128 SSA trajectories, shown by the solid line and shaded area, with the predictions of our momentexpansion method, shown by the dashed line and dotted boundaries.The solid and dashed lines denote average values, while the shaded and dot-bordered areas are the regions within one standard deviation.Panels (a-d) show the statistics of (a) the number of compartments N, (b) the total abundance of the first chemical species M 1,0 , (c) the M 1,1 moment and (d) the Pearson correlation coefficient for the two chemical species.Insets of panels (a,b,c) respectively show the histograms of the N, M 1,0 and M 1,1 moments across the stochastic trajectories, with the black solid line indicating the predicted mean value.The inset of panel (d) shows the 2D histogram of the chemical content of compartments from all trajectories at the final timepoint of the SSA simulation, with the star symbol indicating the Taylor expansion point ( ⟨M 1,0 ⟩ /⟨N⟩, ⟨M 0,1 ⟩ /⟨N⟩).Panels (e-i) show the Pearson correlation and final-timepoint 2D histogram of the system in different k b /k D regimes: from (e) no cell turnover occurring with k b /k D = ∞, increasing the cell turnover speed (f-h) to the fastest (i) with k b /k D = 20.

Remark 1 .
. Each of these products contains, at most, m − 1 moments.Each term in the fully expanded version of eq.(A.3) has the form

Figure 4 :
Figure4: Histograms for the distributions of (a) the cellular content variable x 1 , (b) the number of cells in the system N and (c) the total mass of x 1 in the system M 1 in the Binary birth-death-fusion case study.The solid black lines in each panel shows the estimated values of M 1 / ⟨N⟩, ⟨N⟩ and M 1 , respectively, as computed by the approximation.

2 Figure 6 :
Figure 6: Comparison of the results obtained by averaging 128 SSA trajectories, shown by the solid line and shaded area, with the predictions of our momentexpansion method, shown by the dashed line and dotted boundaries.The solid and dashed lines denote average values, while the shaded and dot-bordered areas are the regions within one standard deviation.The layout of the panels for each row follows that of (fig.3, a-d), while the rows correspond to the different k b /k D regimes as in (fig.3, e-i): (a-d) k b /k D = ∞, (e-h) k b /k D = 20000, (i-l) k b /k D = 2000, (m-p) k b /k D = 200, (q-t) k b /k D = 20.

Table 1 :
Error quantifications in the three case studies: Binary birth-death fusion (bBDF), Shared Antithetic Integral Controller (sAIC) and Mutual gene repression (MR).Errors were calculated by taking the absolute distance between the moments estimated from Monte Carlo simulations and LPAC, dividing it by the Monte Carlo standard deviation and averaging the resulting value across all time points.