Moment-based analysis of biochemical networks in a heterogeneous population of communicating cells

TL;DR

This paper introduces a scalable moment-based method to analyze how cell-cell communication influences noise in biochemical networks within heterogeneous cell populations, balancing accuracy and computational efficiency.

Contribution

It develops a novel, scalable moment-based framework that captures population-level noise effects in communicating cells, using symmetry-based model reduction.

Findings

The method accurately approximates stochastic moments compared to simulations.

Communication can both reduce and amplify noise depending on network parameters.

The approach scales efficiently to large populations without increasing complexity.

Abstract

Cells can utilize chemical communication to exchange information and coordinate their behavior in the presence of noise. Communication can reduce noise to shape a collective response, or amplify noise to generate distinct phenotypic subpopulations. Here we discuss a moment-based approach to study how cell-cell communication affects noise in biochemical networks that arises from both intrinsic and extrinsic sources. We derive a system of approximate differential equations that captures lower-order moments of a population of cells, which communicate by secreting and sensing a diffusing molecule. Since the number of obtained equations grows combinatorially with number of considered cells, we employ a previously proposed model reduction technique, which exploits symmetries in the underlying moment dynamics. Importantly, the number of equations obtained in this way is independent of the number of considered cells such that the method scales to arbitrary population sizes. Based on this approach, we study how cell-cell communication affects population variability in several biochemical networks. Moreover, we analyze the accuracy and computational efficiency of the moment-based approximation by comparing it with moments obtained from stochastic simulations.

Tables

Table 1. Table 1: Normal & lognormal closure functions of order three.

MA	$⟨X_1{X}_2{X}_3⟩$
Normal	$⟨X_1⟩⟨X_2{X}_3⟩+⟨X_2⟩⟨X_1{X}_3⟩$
	$+⟨X_3⟩⟨X_1{X}_2⟩-2⟨X_1⟩⟨X_2⟩⟨X_3⟩$
Lognormal	$\frac{⟨X_1{X}_2⟩⟨X_2{X}_3⟩⟨X_1{X}_3⟩}{⟨X_1⟩⟨X_2⟩⟨X_3⟩}$

Table 2. Table 2: Run times (s) of the reduced moment-based approach and SSA for the positive feedback gene circuit*.

$N$	Reduced moments	SSA (1 sample path)
10	$0.068\pm 0.005$	$0.142\pm 0.004$
100	$0.065\pm 0.004$	$12.300\pm 0.123$
1000	$0.064\pm 0.005$	$2569.20\pm 52.93$

Equations

$$\begin{array}{rl}\sum _{k=1}^S{\alpha }_{j,k}{\mathrm{X}}_{i,k}& \to \sum _{k=1}^S{\beta }_{j,k}{\mathrm{X}}_{i,k}\\ {\mathrm{X}}_{i,S}& \leftrightharpoons {\mathrm{X}}_E,\end{array}$$

$$\begin{array}{rl}\frac{\mathrm{d}P(x,t)}{\mathrm{d}t}& =\sum _{i=1}^N\sum _{j=1}^M[a_{i,j}(x-{\nu }_{i,j},c_{i,j})P(x-{\nu }_{i,j},t)\\ & -a_{i,j}(x,c_{i,j})P(x,t)].\end{array}$$

$$\begin{array}{rl}\frac{\mathrm{d}P(x,t\mid c)}{\mathrm{d}t}& =\sum _{i=1}^N\sum _{j=1}^M[a_{i,j}(x-{\nu }_{i,j},c_{i,j})P(x-{\nu }_{i,j},t\mid c)\\ & -a_{i,j}(x,c_{i,j})P(x,t\mid c)].\end{array}$$

$$\begin{split}\langle\phi(X,C)\rangle&=\big{\langle}\langle\phi(X,C)\mid C\rangle\big{\rangle}\\ &=\Bigg{\langle}\sum_{x\in\mathcal{X}}\phi(x,C)P(x,t\mid C)\Bigg{\rangle},\end{split}$$

$$\small\begin{split}&\frac{\mathrm{d}\langle\phi(X,C)\rangle}{\mathrm{d}t}=\Bigg{\langle}\sum_{x\in\mathcal{X}}\phi(x,C)\sum_{i=1}^{N}\sum_{j=1}^{M}[a_{i,j}(x-\nu_{i,j},C_{i,j})\\ &\qquad P(x-\nu_{i,j},t\mid C)-a_{i,j}(x,C_{i,j})P(x,t\mid C)]\Bigg{\rangle}.\end{split}$$

$$\small\begin{split}&\frac{\mathrm{d}\langle\phi(X,C)\rangle}{\mathrm{d}t}=\Bigg{\langle}\sum_{i=1}^{N}\sum_{j=1}^{M}\Big{\langle}\phi(X+\nu_{i,j},C)a_{i,j}(X,C_{i,j})\mid C\Big{\rangle}\\ &\qquad\qquad\qquad-\Big{\langle}\phi(X,C)a_{i,j}(X,C_{i,j})\mid C\Big{\rangle}\Bigg{\rangle},\end{split}$$

$$\small\begin{split}\frac{\mathrm{d}\langle\phi(X,C)\rangle}{\mathrm{d}t}&=\sum_{i=1}^{N}\sum_{j=1}^{M}\Big{\langle}\phi(X+\nu_{i,j},C)a_{i,j}(X,C_{i,j})\Big{\rangle}\\ &\qquad\qquad\quad-\Big{\langle}\phi(X,C)a_{i,j}(X,C_{i,j})\Big{\rangle}.\end{split}$$

$$K_{eq}=2(NS+1)+\left(\genfrac{}{}{0ex}{}{NS+N{M}^{\prime }+1}{2}\right)-\left(\genfrac{}{}{0ex}{}{N{M}^{\prime }}{2}\right),$$

$$\begin{array}{rl}{\hat{K}}_{eq}=2(S+1)+& \left(\genfrac{}{}{0ex}{}{S+1+M^{\prime }}{2}\right)\\ & +\left(\genfrac{}{}{0ex}{}{S+M^{\prime }}{2}\right)-2\left(\genfrac{}{}{0ex}{}{M^{\prime }}{2}\right)+S,\end{array}$$

$$\begin{array}{c}\begin{array}{r}{\mathrm{X}}_{i,1}\stackrel{C_{i,1}}{\to }{\mathrm{X}}_{i,2}\end{array}\begin{array}{r}{\mathrm{X}}_{i,2}\underset{c_3}{\overset{c_2}{\leftrightharpoons }}{\mathrm{X}}_E,\end{array}\end{array}$$

$$\frac{\mathrm{d}⟨X_{i,2}⟩}{\mathrm{d}t}=⟨C_{i,1}{X}_{i,1}⟩+c_2⟨X_E⟩-c_3⟨X_{i,2}⟩.$$

$$\begin{array}{rl}\frac{\mathrm{d}⟨X_{i,2}{X}_E⟩}{\mathrm{d}t}& =⟨C_{i,1}{X}_{i,1}{X}_E⟩+c_2⟨X_E^2⟩-c_2⟨X_E⟩\\ & +c_3⟨X_{i,2}^2⟩-c_3⟨X_{i,2}{X}_E⟩-c_3⟨X_{i,2}⟩\\ & +\sum _{k\ne i}{c}_3⟨X_{k,2}{X}_{i,2}⟩-c_{2}N⟨X_{i,2}{X}_E⟩.\end{array}$$

$$\begin{array}{rl}\frac{\mathrm{d}⟨X_{i,2}{X}_E⟩}{\mathrm{d}t}& =⟨C_{i,1}{X}_{i,1}{X}_E⟩+c_2⟨X_E^2⟩-c_2⟨X_E⟩\\ & +c_3⟨X_{i,2}^2⟩-c_3⟨X_{i,2}{X}_E⟩-c_3⟨X_{i,2}⟩\\ & +c_3(N-1)⟨X_{j,2}{X}_{i,2}⟩-c_{2}N⟨X_{i,2}{X}_E⟩.\end{array}$$

$$\begin{array}{c}\begin{array}{r}\varnothing \stackrel{C_i^b}{\to }{\mathrm{P}}_i\stackrel{C_i^d}{\to }\varnothing \end{array}\begin{array}{r}{\mathrm{P}}_i\underset{c^t}{\overset{c^t}{\leftrightharpoons }}\mathrm{Q},\end{array}\end{array}$$

$$\mathrm{CV}[P_i]=\sqrt{\frac{⟨P_i^2⟩-⟨P_i{⟩}^2}{⟨P_i{⟩}^2}},$$

$$\mathrm{PV}[P_i,P_k]=\sqrt{\frac{⟨(P_i-P_k{)}^2⟩}{⟨P_i⟩⟨P_k⟩}},$$

$$\begin{array}{c}\begin{array}{r}\varnothing \stackrel{C_i^b}{\to }{\mathrm{A}}_i\stackrel{C_i^d}{\to }\varnothing \\ {\mathrm{A}}_i\stackrel{c^a}{\to }{\mathrm{A}}_i+{\mathrm{A}}_i\end{array}\begin{array}{r}{\mathrm{A}}_i\underset{c^t}{\overset{c^t}{\leftrightharpoons }}\mathrm{B}.\end{array}\end{array}$$

$$\begin{array}{c}\begin{array}{r}{\mathrm{D}}_i\stackrel{C_i^m}{\to }{\mathrm{D}}_i+{\mathrm{M}}_i\\ {\mathrm{M}}_i\stackrel{C_i^p}{\to }{\mathrm{M}}_i+{\mathrm{P}}_i\\ {\mathrm{D}}_i+{\mathrm{P}}_i\underset{c^d}{\overset{c^a}{\leftrightharpoons }}{\mathrm{DP}}_i\\ {\mathrm{DP}}_i\stackrel{c^f}{\to }{\mathrm{DP}}_i+{\mathrm{M}}_i\end{array}\begin{array}{r}{\mathrm{M}}_i\stackrel{d^m}{\to }\varnothing \\ {\mathrm{P}}_i\stackrel{d^p}{\to }\varnothing \\ {\mathrm{P}}_i\underset{c^t}{\overset{c^t}{\leftrightharpoons }}\mathrm{Q},\end{array}\end{array}$$

Full text

Moment-based analysis of biochemical networks in a heterogeneous population of communicating cells

David T. Gonzales

Max Planck Institute of Molecular Cell Biology and Genetics, 01307 Dresden, Germany

Center for Systems Biology Dresden, 01307 Dresden, Germany

T-Y Dora Tang

Max Planck Institute of Molecular Cell Biology and Genetics, 01307 Dresden, Germany

Christoph Zechner

Max Planck Institute of Molecular Cell Biology and Genetics, 01307 Dresden, Germany

Center for Systems Biology Dresden, 01307 Dresden, Germany

Correspondence to: zechner@mpi-cbg.de

**Cells can utilize chemical communication to exchange information and coordinate their behavior in the presence of noise. Communication can reduce noise to shape a collective response, or amplify noise to generate distinct phenotypic subpopulations. Here we discuss a moment-based approach to study how cell-cell communication affects noise in biochemical networks that arises from both intrinsic and extrinsic sources. We derive a system of approximate differential equations that captures lower-order moments of a population of cells, which communicate by secreting and sensing a diffusing molecule. Since the number of obtained equations grows combinatorially with number of considered cells, we employ a previously proposed model reduction technique, which exploits symmetries in the underlying moment dynamics. Importantly, the number of equations obtained in this way is independent of the number of considered cells such that the method scales to arbitrary population sizes. Based on this approach, we study how cell-cell communication affects population variability in several biochemical networks. Moreover, we analyze the accuracy and computational efficiency of the moment-based approximation by comparing it with moments obtained from stochastic simulations. **

1  Introduction

In recent years, significant progress has been made in developing computational models and algorithms to study stochastic biochemical networks inside living cells. Most commonly, these models are based on the chemical master equation (CME), whose solution provides a time-dependent probability distribution over molecular concentrations [1]. CME-based models can account for the discrete and random nature of biochemical reactions (intrinsic noise) as well as population heterogeneity stemming from differences in each cell’s microenvironment (extrinsic variability) [2]. The computational analysis of the CME is challenging but several efficient numerical techniques have been proposed in the past, including stochastic simulation algorithms [3], moment-based methods [4, 5] and combinations thereof [6, 7, 8].

However, the majority of existing approaches to study noise in cell populations rely on the assumption that individual cells act independently of each other. More concretely, each cell’s dynamics is considered to be an independent and identically distributed realization of the same stochastic process. Evidently, this assumption is violated in systems where cells communicate with one another to coordinate their behavior [9, 10]. Typical examples include quorum-sensing systems in bacterial colonies [11], or paracrine communication in higher organisms [12]. Understanding the interplay between cell-cell communication and the stochastic behavior of individual cells is an important challenge and demands for suitable mathematical approaches. However, extending standard techniques to account for cell-cell communication leads to computational difficulties, because the dimensionality of the resulting models increases with the number of cells in a population.

Recently, first attempts have been made to develop more tractable stochastic models of systems of communicating cells. In [9], for instance, the authors use a moment-based method to study how neighbour-neighbour-coupling affects concentration fluctuations in a tissue. A related approach has been proposed in [13] to study community effects in cells that interact with each other by secreting and sensing certain signalling molecules. In particular, the authors show how the dimensionality of the considered model can be dramatically reduced by exploiting certain symmetries in the governing equations. However, both approaches account exclusively for intrinsic noise, whereas extrinsic sources have not been considered. In [14], the authors study how chemical communication via a quorum sensing molecule affects intrinsic and extrinsic noise in cell communities. To obtain tractable simulations, they used a tailored approach that combines stochastic simulations with a quasi-steady state approximation to eliminate fast variables from the model.

In this article, we discuss a moment-based approach to study noise in communicating cells that arises from both intrinsic and extrinsic sources. In particular, we focus on a secrete-and-sense model [10, 13], where individual cells can secrete a signalling molecule to the extracellular environment, which in turn can be sensed by other cells in the population. We first derive a system of differential equations that captures lower order moments of the population. Since the number of obtained equations grows combinatorially with number of considered cells, we perform a symmetry-based model reduction based on the work of Batmanov et al. (2012) [13]. The dimensionality of the reduced model is independent of the population size, which makes it computationally very efficient. We employ this approach to study how noise is affected by cell-to-cell communication in a few biochemical systems.

The rest of the paper is structured as follows. A general model of biochemical networks in chemically interacting cell populations is presented in Section 2. In Section 3, we provide a stochastic description of this model based on the CME and derive a general equation for its moment dynamics. In Section LABEL:reduction, we show how the number of moment equations can be reduced by exploiting symmetries of the model. Finally, we apply our approach to analyze several biochemical networks in Section 5.

2  Reaction networks in a population of chemically interacting cells

We consider a population of $N$ genetically identical cells that communicate with each other through a diffusing signalling molecule (Fig. 1). Each individual cell $i$ is associated with an identical set of $S$ chemical species $\mathrm{X}_{i,1},\ldots,\mathrm{X}_{i,S}$ that interact with one another via $M$ biochemical reaction channels. Without loss of generality, we consider the first $S-1$ species to be confined to the intracellular environment of cell $i$. The $S$th species corresponds to the signalling molecule, which can shuttle between the intra- and extracellular environment through transport reactions. For simplicity, we consider the external environment to be well-mixed, such that the import of signalling molecules into any cell $i$ does not depend on the spatial configuration of the system. In total, the system can be described by a reaction network

[TABLE]

for $i=1,\ldots,N$ and $j=1,\ldots,M-2$. In (1), $\alpha_{j,k}$ and $\beta_{j,k}$ correspond to the reactant and product coefficients of the respective reaction. The species $\mathrm{X}_{E}$ denotes the signalling molecule in the external environment.

In total, the network comprises $NS+1$ chemical species and $NM$ chemical reactions. We define by $X_{i}(t)=(X_{i,1}(t),\ldots,X_{i,S}(t))$ the state of cell $i$, which collects the copy numbers of all species associated with this cell at time $t$. The state of the overall system is then given by $X(t)=(X_{1}(t),\ldots,X_{N}(t),X_{E}(t))$, with $X_{E}(t)$ as the number of signalling molecules in the external environment. Moreover, we denote by $\nu_{i,k}\in\mathbb{Z}^{NS+1}$ the stoichiometric change vector associated with the $k$th reaction channel of cell $i$ consistent with reaction network (1).

3  Moment dynamics of heterogeneous cell communities

If both the intra- and extracellular environments are well-mixed, we can describe the state $X(t)$ as a continuous-time Markov chain, whose state probability distribution $P(x,t)=P(X(t)=x)$ admits a master equation of the form

[TABLE]

In (2), the function $a_{i,j}$ is the reaction propensity associated with reaction $j$ in cell $i$. Throughout this article, we consider the propensities to obey the law of mass action such that $a_{i,j}(x,c_{i,j})=c_{i,j}g_{j}(x_{i},x_{E})$ with $x_{i}$ as the part of the state vector associated with cell $i$, $x_{E}$ as the abundance of the signalling molecule in the environment, $c_{i,j}\in\mathbb{R}^{+}$ as a stochastic rate constant, and $g_{j}$ as a polynomial. Note that while $g_{j}$ is identical for all cells, we allow the individual rate constants $c_{i,j}$ to vary across the population. This provides a means to account for extrinsic sources of cell-to-cell variability [5]. We consider the reaction rate constants for each cell $i$ to be independent random vectors $C_{i}=(C_{i,1},\ldots,C_{i,M})$ drawn from a common probability distribution $C_{i}\sim p_{c}(\cdot)$ for $i=1,\ldots,N$. Note that deterministic (non-varying) reaction rates can be accounted for by letting $p_{c}$ be a Dirac measure with respect to this parameter. With $C=(C_{1},\ldots,C_{N})$, we can formulate a master equation for the conditional distribution $P(x,t\mid c)=P(X(t)=x\mid C=c)$, i.e.,

[TABLE]

3.1  Moment dynamics

To analyze the cell community model, we resort to a moment-based approach, which provides a lower-dimensional description of the population and its heterogeneity. More precisely, we seek for the population moments

[TABLE]

with $\phi:(x,c)\rightarrow\mathbb{R}$ as a monomial in $(x,c)$ and $\mathcal{X}$ as the domain of $X(t)$. Note that we omit the dependency of the moments on time for the sake of a compact notation. In order to derive a differential equation for the time evolution of $\langle\phi(X,C)\rangle$, we calculate the derivative of (4) and insert the r.h.s. of (3)

[TABLE]

Using a change of variable, it is straightforward to show that (5) simplifies to

[TABLE]

where the inner brackets denote expectations conditionally on the random parameters $C$. Now, using double expectations, we obtain

[TABLE]

Eq. (7) describes the time evolution of moments and cross moments for a heterogeneous population of secrete-and-sense cells. It can thus be seen as an extension of the moment equations derived in [5] to account for cell-to-cell communication. Indeed, if we set the transport rates to zero, all cells in the population become independent of each other such that the two approaches yield equivalent solutions.

Note that depending on the details of the system, the moment equations from (7) may not be closed. This is the case in the presence of second-order reactions, or even first-order reactions if their rate constants are randomly distributed across the population. Moment-closure approximation (MA) techniques provide a popular means to address this problem by imposing certain assumptions on the underlying state probability distribution [4, 5, 15]. More precisely, this allows us to replace the higher order moments that appear on the r.h.s. of (7) by functions of the lower order moments. These functions are referred to as closure functions and their particular form depends on the distributional assumption we make. Popular choices include the normal [15] and lognormal [4] closure functions and we will adopt those in the present study. To check the accuracy of our MAs in this study, we compare them with moments calculated using the Stochastic Simulation Algorithm (SSA) [3]. Throughout this article, we consider moments of up to second order and replace all third-order moments that appear on the r.h.s. of (7) using MA functions provided in Table 1.

4  Symmetry-based model reduction

While eq. (7) provides a more tractable description of the cell community than (2), its dimensionality still scales combinatorially with the number of considered cells $N$. In case we consider all first and second order moments, the number of equations is given by

[TABLE]

where $M^{\prime}$ is the number of rate constants among the $M$ reactions that vary from cell to cell. The first term in (8) accounts for 1st and 2nd order moments of the chemical species, and the remaining terms account for all cross-moments of species and rate constants that change over time.

The number of first and second order moment equations scales quadratically with $N$, which limits the above approach to relatively small population sizes. However, the moment system can be reduced substantially by taking into account the symmetries of the considered model as has been proposed in [13]. More precisely, if we consider all initial cell states $X_{i}(0)$ to be identically distributed, the moment dynamics of each cell will be equivalent and indistinguishable for all times $t>0$ such that $\langle X_{i,k}\rangle=\langle X_{j,k}\rangle$, $\langle X_{i,k}^{2}\rangle=\langle X_{j,k}^{2}\rangle$, $\langle X_{i,k}X_{j,l}\rangle=\langle X_{m,k}X_{n,l}\rangle$, and $\langle X_{i,k}X_{E}\rangle=\langle X_{j,k}X_{E}\rangle$ for any $i\neq j$, $m\neq n$, $k$ and $l$. Consequently, we can obtain a reduced model by considering only the moments up to the second order and cross-moments of the states of any two reference cells $X_{i}(t)$ and $X_{j}(t)$ as well as the amount of signalling molecules in the external environment $X_{E}(t)$. The resulting set of equations can be reduced further by eliminating all moments associated with one of the two cells $j$ (e.g., $\langle X_{j,k}\rangle$ or $\langle X_{j,k}X_{E}\rangle$) since those are identical to the corresponding moments of cell $i$. In total, the required number of equations is then given by

[TABLE]

which is thus independent of the population size.

4.1  Illustrative example

Consider the toy model

[TABLE]

with $i=1,\ldots,N$. We consider the rate $C_{i,1}$ to be randomly distributed across the population, while $c_{2}$ and $c_{3}$ are identical in all cells.

To demonstrate how the original system can be reduced based on symmetries, we distinguish between two different cases. The first case concerns equations that are the same in the original and reduced model. These are equations of the moments that do not involve the signalling molecule in the external environment such as $\langle X_{i,1}\rangle$, $\langle X_{i,2}\rangle$, or $\langle X_{i,1}X_{i,2}\rangle$. For example, in both the original and reduced models, the expectation of species $X_{i,2}$ satisfies

[TABLE]

The dynamics of moments and cross-moments involving the external signalling molecule depend on all cells in the population due to the transport reactions. In the case of $\langle X_{i,2}X_{E}\rangle$, the original equation is:

[TABLE]

The terms in the sum capture the dependencies between two cells $k$ and $i$. Since all terms in the sum are identical due to the symmetry of the population, they can be replaced by the contribution of any cell $j\neq i$ multiplied by $(N-1)$,

[TABLE]

Therefore, if we perform analogous manipulations for all other moments and cross-moments involving the signalling molecule (e.g., $\langle X_{E}\rangle$, $\langle X_{i,1}X_{E}\rangle$) and eliminate all remaining redundancies between cell $i$ and $j$, we arrive at a system of $17$ coupled differential equations, independent of the population size $N$. For further information on the symmetry-based model reduction approach, the reader may refer to [13].

5  Case studies

In this section, we use the described moment-based approach to study how cell-cell communication affects noise in different biochemical networks. Python scripts used for this study are available in github.com/zechnerlab/CommunityMoments. Moment equations were solved using a SciPy numerical solver (solve_ivp). SSA sample paths from Gillespie’s Direct Method were simulated using Tellurium [16] and used to test the accuracy of the moment-based approach. For the sake of a compact notation, molecular species are assigned different letters and reaction rate parameters are assigned letters with superscripts.

5.1  Birth-death process

As a first example, we study a birth-death process in a heterogeneous population of interacting cells, i.e.,

[TABLE]

for $i=1,\ldots,N$. Here, the birth and death reaction rate constants are independent random variables as indicated by capital letters $C^{b}_{i}$ and $C^{d}_{i}$, while the transport rate $c^{t}$ is fixed. For this reaction network, a symmetry-reduced model of up to the second order can be described using $\hat{K}_{eq}=12$ differential equations regardless of the population size.

The goal of this case study is to study how cell-to-cell communication affects the variability in the abundance of $\mathrm{P}_{i}$. To quantify variability, we use two metrics: the first one is the coefficient of variation (CV) defined as

[TABLE]

for any $i=1,\ldots,N$. The CV captures the expected variation in protein abundance inside single cells across different populations. Due to the symmetry, we have $\mathrm{CV}[P_{i}]=\mathrm{CV}[P_{k}]$ for any $i$ and $k$. Furthermore, we define the pair variation (PV)

[TABLE]

which captures the expected variation between two different cells $i$ and $k$ within the same population. Note that in the absence of cell-cell communication, (15) and (16) are identical up to a scaling factor of $\sqrt{2}$. For a communicating population, however, this is not the case due to correlations between cells in the population.

In Figs. 2A and B, we compare the first and second order moments of species $P_{i}$ for the symmetry-reduced model of size $N=10$ with stochastic simulations and found a good agreement between the two approaches. In Figs. 2C-F, we show the dependency of CV[$P_{i}$] and PV[$P_{i}$,$P_{j}$] as a function of the transport rate as well as the population size. The coefficient of variation CV[$P_{i}$] decreases with increasing transport rates and to some extent also with the population size. While the pair variation PV[$P_{i}$,$P_{j}$] shows a similar decrease with increasing transport rates, it seems to be independent of the population size.

We also show how extrinsic variability changes the steady-state variability of $\mathrm{P}_{i}$ by increasing the CV of the birth rate parameter $C^{b}_{i}$, while keeping all other parameters constant. Both CV[$P_{i}$] and PV[$P_{i}$,$P_{j}$] increase with increasing extrinsic variability, whereas large transport rates can attenuate this effect due to spatial averaging of $P_{i}$ levels (Fig. 3).

5.2  Autocatalytic circuit

Next, we focus on an autocatalytic system defined by

[TABLE]

Here we consider $C^{b}_{i}$ and $C^{d}_{i}$ to be randomly distributed across the population, while the autocatalytic rate $c^{a}$ and transport rate $c^{t}$ are fixed. In this example, the total number of reduced equations is $\hat{K}_{eq}=12$. To obtain a closed set of moments, we applied the lognormal closure and checked its accuracy by comparing it to Monte Carlo estimates of the moments calculated over $1000$ SSA realizations (Fig. 4). We generally found a good agreement between the MA and the SSA. Similar to the birth-death systems, we observe that both CV and PV of species $\mathrm{A}_{i}$ decrease with increasing transport rates. Increasing the population size does not affect the PV, but decreases its CV.

To study the relationship between the CV and PV, we plotted CV[$A_{i}$] against PV[$A_{i}$,$A_{j}$] in Figure 5. As expected, CV and PV are related by a factor of $\sqrt{2}$ when $c^{t}=0$. Slight deviations from this scaling are possible due to the approximations involved in the derivation of the moment equations. In the presence of communication, PV[$A_{i}$,$A_{j}$] drops below the $\sqrt{2}$ scaling law, indicating that the variability between cells in the same population is smaller than the variability of cells across different populations.

5.3  Genetic feedback circuit

Lastly, we tested the moment-based method using a larger system. In particular, we focus on a genetic feedback circuit given by

[TABLE]

for $i=1,\ldots,N$. Here, we consider the reaction rate constants associated with transcription $C^{m}_{i}$ and translation $C^{p}_{i}$ to be randomly distributed across the population, while all other rate parameters are fixed. Transcriptional feedback of the gene circuit is mediated by the protein product ($\mathrm{P}$), which binds DNA ($\mathrm{D}$) to form a complex $\mathrm{DP}$. In the bound state, the gene can be transcribed with rate constant $c^{f}$. Depending on the ratio of the bound and unbound transcription rate, the feedback mechanism can either enhance ($c^{f}>\langle C^{m}_{i}\rangle$) or inhibit ($c^{f}<\langle C^{m}_{i}\rangle$) gene expression. We consider a population of $N=10$ cells and applied a lognormal closure to solve for the reduced moment dynamics. The total number of reduced equations for this system is $\hat{K}_{eq}=48$. Fig. 6 shows that the time evolution of the mRNA moments obtained from the symmetry-reduced model and SSA are very similar to each other.

To analyze the computational efficiency of the moment-based approach, we recorded run times for the positive feedback gene circuit for different population sizes $N$ and compared it against the SSA (Table 2). As expected, the run time for the moment-based model does not increase with $N$ while the run time associated with the SSA increases by orders of magnitude.

6  Conclusions

We presented and analyzed an efficient moment-based approach to study noise in heterogeneous populations of communicating cells. By exploiting certain symmetries of the underlying moment dynamics [13], the number of differential equations needed to describe the population could be strongly reduced. Indeed, the dimensionality of the reduced model is independent of the number of considered cells, which enables the analysis of regulatory processes in large populations. In principle, the method applies to arbitrary intracellular reaction networks, although its accuracy relies on the availability of suitable moment closure functions. In the present work, we tested the approach using three different case studies, where it accurately captured first and second order moments using normal and lognormal closure functions. Our analyses show that moment-based approximations provide an effective way to study how biochemical noise and population heterogeneity affect the collective dynamics of communicating cells.