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

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.

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

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

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

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

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

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

which is thus independent of the population size.

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.

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.