The Community Simulator: A Python package for microbial ecology

TL;DR

The Community Simulator is an open-source Python package that enables scalable, reproducible simulations of microbial community dynamics, incorporating advanced algorithms for equilibrium finding and supporting complex ecological modeling.

Contribution

It introduces a comprehensive Python toolkit with novel EM algorithms for efficient equilibrium computation and flexible modeling features for microbial ecology simulations.

Findings

EM algorithm improves performance by 10-100x over direct integration

Supports large-scale, reproducible microbial community simulations

Facilitates analysis of microbiome data and ecological dynamics

Abstract

Natural microbial communities contain hundreds to thousands of interacting species. For this reason, computational simulations are playing an increasingly important role in microbial ecology. In this manuscript, we present a new open-source, freely available Python package called Community Simulator for simulating microbial population dynamics in a reproducible, transparent and scalable way. The Community Simulator includes five major elements: tools for preparing the initial states and environmental conditions for a set of samples, automatic generation of dynamical equations based on a dictionary of modeling assumptions, random parameter sampling with tunable levels of metabolic and taxonomic structure, parallel integration of the dynamical equations, and support for metacommunity dynamics with migration between samples. To significantly speed up simulations using Community Simulator, our Python package implements a new Expectation-Maximization (EM) algorithm for finding equilibrium states of community dynamics that exploits a recently discovered duality between ecological dynamics and convex optimization. We present data showing that this EM algorithm improves performance by between one and two orders compared to direct numerical integration of the corresponding ordinary differential equations. We conclude by listing several recent applications of the Community Simulator to problems in microbial ecology, and discussing possible extensions of the package for directly analyzing microbiome compositional data.

Tables

Table 1. Table 1: Parameters and units for the Microbial Consumer Resource Model.

$N_i$	population density of species $i$ (individuals/volume)
$R_{\alpha }$	Concentration of resource $\alpha$ (mass/volume)
$c_{i\alpha }$	Uptake rate per unit concentration of resource $\alpha$ by species $i$ (volume/time)
$D_{\alpha \beta }$	Fraction of byproducts from resource $\beta$ converted to $\alpha$ (unitless)
$g_i$	Conversion factor from energy uptake to growth rate (1/energy)
$w_{\alpha }$	Energy content of resource $\alpha$ (energy/mass)
$l_{\alpha }$	Leakage fraction for resource $\alpha$ (unitless)
$m_i$	Minimal energy uptake for maintenance of species $i$ (energy/time)
$R_{\alpha }^0$	Intrinsic equilibrium abundance of resource $\alpha$ (mass/volume)
${\tau }_{\alpha }$	Timescale for externally supplied resource turnover (time)
$r_{\alpha }$	Rate of resource self-renewal (volume/mass/time)
$n$	Hill coefficient for functional response (unitless)
${\sigma }_{\max }$	Maximum input flux (mass/time)
$n_{\mathrm{reg}}$	Hill coefficient for metabolic regulation (unitless)

Table 2. Table 2: Definitions of global parameters used for constructing random ecosystems. Values of these parameters are supplied as a Python dictionary to the function MakeMatrices , which generates randomly sampled consumer preference and metabolic matrices.

$M$	Number of resources
$T$	Number of resource classes
$S_{\mathrm{tot}}$	Number of microbial species in regional pool
$F$	Number of specialist families
$S$	Number of microbial species initially present in each local community
${\mu }_c$	Mean sum over a row of the preference matrix $c_{i\alpha }$
${\sigma }_c$	Standard deviation of sum over a row for Gaussian or Gamma $c_{i\alpha }$
$c_0$	Low consumption level for Binary $c_{i\alpha }$
$c_1$	High consumption level for Binary $c_{i\alpha }$
$q$	Fraction of consumption capacity allocated to preferred resource class
$s$	Sparsity of metabolic matrix
$f_w$	Fraction of secreted byproducts allocated to “waste” resource class
$f_s$	Fraction of secreted byproducts allocated to same resource class

Equations

$$J^{\mathrm{in}}=J^{\mathrm{growth}}+J^{\mathrm{out}}$$

$${\nu }_{\beta }^{\mathrm{out}}\equiv J_{\beta }^{\mathrm{out}}/w_{\beta }$$

$${\nu }_{\beta }^{\mathrm{in}}\equiv J_{\beta }^{\mathrm{in}}/w_{\beta }$$

$$\frac{d{N}_i}{dt}=g_i{N}_i(J_i^{\mathrm{growth}}-m_i).$$

$$\frac{d{R}_{\alpha }}{dt}=h_{\alpha }(R_{\alpha })-\sum _j{N}_j{\nu }_{j\alpha }^{\mathrm{in}}+\sum _j{N}_j{\nu }_{j\alpha }^{\mathrm{out}},$$

$$h_{\alpha }^{\mathrm{external}}(R_{\alpha })={\tau }_{\alpha }^{-1}(R_{\alpha }^0-R_{\alpha })$$

$$h_{\alpha }^{\mathrm{self}-\mathrm{renewing}}(R_{\alpha })=r_{\alpha }{R}_{\alpha }(R_{\alpha }^0-R_{\alpha }).$$

$${\nu }_{i\alpha }^{\mathrm{in}}=\sigma (c_{i\alpha }{R}_{\alpha })$$

$${\sigma }_I(x)=x,$$

$${\sigma }_{II}(x)=\frac{x}{1+\frac{x}{{\sigma }_{\max }}}$$

$${\sigma }_{III}(x)=\frac{x^n}{1+\frac{x^n}{{\sigma }_{\max }}},$$

$$J^{\mathrm{out}}=l{J}^{\mathrm{in}}.$$

$$J_i^{\mathrm{growth}}=\sum _{\alpha }(1-l_{\alpha })J_{i\alpha }^{\mathrm{in}}=\sum _{\alpha }(1-l_{\alpha })w_{\alpha }\sigma (c_{i\alpha }{R}_{\alpha })$$

$$J_{i\beta }^{\mathrm{out}}=\sum _{\alpha }{D}_{\beta \alpha }{l}_{\alpha }{J}_{i\alpha }^{\mathrm{in}}=\sum _{\alpha }{D}_{\beta \alpha }{l}_{\alpha }{w}_{\alpha }\sigma (c_{i\alpha }{R}_{\alpha }).$$

$${\nu }_{i\beta }^{\mathrm{out}}=\sum _{\alpha }{D}_{\beta \alpha }{l}_{\alpha }\frac{w_{\alpha }}{w_{\beta }}\sigma (c_{i\alpha }{R}_{\alpha })$$

$$\frac{d{N}_i}{dt}$$

$$\frac{d{R}_{\alpha }}{dt}$$

$$u_{i\alpha }^{\mathrm{in}-w}=\frac{(w_{\alpha }{c}_{i\alpha }{R}_{\alpha }{)}^{n_{\mathrm{reg}}}}{{\sum }_{\beta }(w_{\beta }{c}_{i\beta }{R}_{\beta }{)}^{n_{\mathrm{reg}}}},$$

$$u_{i\alpha }^{\mathrm{in}-\nu }=\frac{(c_{i\alpha }{R}_{\alpha }{)}^{n_{\mathrm{reg}}}}{{\sum }_{\beta }(c_{i\beta }{R}_{\beta }{)}^{n_{\mathrm{reg}}}}$$

$${\nu }_{\beta }^{\mathrm{in}}=u_{i\beta }^{\mathrm{in}}\sigma (c_{i\beta }{R}_{\beta })$$

$$\frac{d{N}_i}{dt}$$

$$\frac{d{R}_{\alpha }}{dt}$$

$$⟨c_{i\alpha }^{\mathrm{gen}}⟩$$

$$⟨(\delta c_{i\alpha }^{\mathrm{gen}}{)}^2⟩$$

$$⟨c_{i\alpha }^A⟩=\{\begin{array}{ll}\frac{{\mu }_c}{M}\left[1+\frac{M-M_A}{M_A}q\right],& \text{if}\alpha \in \mathbf{A}\\ \frac{{\mu }_c}{M}(1-q),& \text{otherwise},\end{array}$$

$$⟨(\delta c_{i\alpha }^A{)}^2⟩=\{\begin{array}{ll}\frac{{\sigma }_c^2}{M}\left[1+\frac{M-M_A}{M_A}q\right],& \text{if}\alpha \in \mathbf{A}\\ \frac{{\sigma }_c^2}{M}(1-q),& \text{otherwise}.\end{array}$$

$$c_{i\alpha }^A=\frac{c_0}{M}+c_1{X}_{i\alpha },$$

$$p_{i\alpha }^A=\{\begin{array}{ll}\frac{{\mu }_c}{M{c}_1}\left[1+\frac{M-M_A}{M_A}q\right],& \text{if}\alpha \in A\\ \frac{{\mu }_c}{M{c}_1}(1-q),& \text{otherwise}\end{array}$$

$$p_{i\alpha }^{\mathrm{gen}}=\frac{{\mu }_c}{M{c}_1}$$

$$D_{\alpha \beta }$$

$$d_{\alpha \beta }$$

$$d^{\mathrm{external}}({\mathbf{R}}^0,\mathbf{R})=\sum _{\beta }{w}_{\beta }{\tau }_{\beta }^{-1}\left[R_{\beta }^0\ln \frac{R_{\beta }^0}{R_{\beta }}-(R_{\beta }^0-R_{\beta })\right]$$

Full text

Background

The last decade has seen a renewed interest in the study of microbial communities. Different environments can harbor diverse communities containing from hundreds to thousands of distinct microbes EMP ; HMP . A central goal of community ecology is to understand the ecological processes that shape these diverse ecosystems. The diversity and function of ecosystems are affected by a wide variety of factors including energy and resource availability Loreau1995 ; Embree2015 , ecological processes such as competition between speciesGause1935 ; MacArthur1970 ; Levin1970 ; Chesson1990 and stochastic colonization Chase2003 ; Jeraldo2012 ; Kessler2015 ; Vega2017 .

Microbial ecosystems also present several new challenges specific to microbes that are not usually addressed in the theoretical ecology literature. Classical models of community ecology (especially niche-based theories) have traditionally considered ecosystems with a few species and resources tilman_resource_1982 ; chesson2000mechanisms . However, microbial ecosystems often have thousands of species and hundreds of small molecules that can be consumed. It is unclear how the intuitions and results from these low-dimensional settings scale to microbiomes. It is known that diverse ecosystems can exhibit distinct emergent features and phase transitions not found in low-dimensional systems Fisher2014 ; Dickens2016 ; bunin2017ecological ; Barbier2018 . Furthermore, classical ecological models usually assume a strict trophic layer separation, ignoring cross-feeding and syntrophy – the consumption of metabolic byproducts of one species by another species. It is now becoming clear that cross-feeding is a central component of microbial ecosystems pacheco2019costless ; Goldford2018 ; muscarella2019species ; marsland2018available and any ecological model must account for this phenomenon.

For these reasons, there is a need for new ways of understanding microbial ecosystems. One powerful approach for understanding complex systems is through simulations. However, simulating diverse microbial ecosystems presents some unique challenges. First, most ecosystems are mathematically represented by complicated coupled, non-linear ordinary differential equations. Simulating these systems in ecosystems with hundreds to thousands of species and metabolites becomes computationally difficult and time-consuming. Second, these dynamical models have thousands of parameters. One needs a principled and biologically realistic way of choosing such parameters. Third, explaining real data requires incorporating ecological processes such as stochastic colonization that play an important role in shaping community structure and dynamics. Finally, we need to be able to incorporate spatial and population level structures in an experimentally realistic way.

Recently, we presented a powerful minimal model of microbial ecosystems that addresses these concerns Goldford2018 ; marsland2018available . Furthermore, we have found a mathematical mapping between ecological dynamics and constrained optimization that can be used to accelerate simulations of many large ecosystems mehta2018constrained ; Marsland2019a . In this paper, we present a new open-source Python package for microbial ecology called Community Simulator that implements these theoretical advances, making it easy to simulate complex microbial communities in a variety of experimentally relevant settings.

Implementation

The architecture of Community Simulator is inspired by the parallel experiments commonly performed with 96-well plates, as illustrated in Fig 1. The central object of the package is a Community class, whose instances are initialized by specifying the initial population sizes and resource concentrations for each parallel “well,” along with the functions and parameters that define the population dynamics. The initial state, dynamical equations and parameters can all be generated automatically from a dictionary of modeling assumptions, or custom-built by the user. Each instance of this class represents an $n$-well plate, containing $n$ well-mixed, non-interacting communities. Once initialized, the state of the plate can be updated in one of two ways. Propagate(T) propagates the system for a time $T$ by integrating the supplied dynamical equations, and Passage(f) builds a replacement plate by adding a fraction $f_{\mu\nu}$ of the contents of each old well $\nu$ to each new well $\mu$. In the final section below, we will discuss a third method SteadyState(), which can find the fixed point of the dynamics in some models without numerical integration.

The package also includes some functions for analysis of the simulation results, including a variety of measures of alpha diversity, as well as extraction of energy flux networks, effective interaction coefficients, and sensitivities to parameter perturbations.

In the following sections, we describe the functionality of each element of the package in turn. We particularly focus on the tools for generating the dynamical equations and the parameter sets, explaining how increasing levels of biological realism can be progressively incorporated. Each section has a corresponding segment in the Jupyter notebook Tutorial.ipynb included with the Community Simulator package. This notebook contains all the code and parameters for generating the figures found in the paper.

Constructing the initial state

The state of a Community instance is contained in a pair of Pandas data frames (https://pandas.pydata.org/, mckinney-proc-scipy-2010 ), one of size $S_{\rm tot}\times n$ for the microbial population sizes $N_{i}$ ($i=1,2,\dots S_{\rm tot}$) and one of size $M\times n$ for the resource abundances $R_{\alpha}$ ($\alpha=1,2,\dots M$). Each row of the data frame corresponds to a different species or resource type, while each column corresponds to a different well.

The function MakeInitialState automatically creates data frames N0,R0 of initial population sizes and resource abundances corresponding to some common experimental scenarios, specified in a dictionary of assumptions. The initial species abundances are supplied by a stochastic process that is agnostic to species identity. This roughly captures the various dispersal mechanisms including mechanical disturbances and turbulent flow that convey microbial cells to new environments. Specifically, random subsets of $S$ species from the regional pool of size $S_{\rm tot}$ are supplied to the $n$ wells of the plate. The population sizes of these species are set to 1 by default, and can be rescaled afterwards if desired.

Initial resource abundances are generated by MakeInitialState based on a Biolog plate scenario, where each well is supplied with a single carbon source. The assumptions dictionary specifies the identity and quantity of the carbon source for each well. Arbitrarily chosen resource abundances can of course be directly supplied to the Community instance instead, to simulate more general conditions.

To capture coarse-grained metabolic structure, the $M$ resources can be assigned to $T$ classes (e.g. sugars, amino acids, etc.), each with $M_{A}$ resources where $A=1,\ldots T$ and $\sum_{A}M_{A}=M$. These class labels become functionally relevant by biasing the sampling of consumption preferences and byproduct stoichiometry, as will be described below. Likewise the $S_{\rm tot}$ species can be assigned to $F$ families, with $F\leq T$, and each family preferentially consuming resources from a different resource class. A generalist family can also be included, with $S_{\rm gen}$ species and no preferred resource class, so that $S_{\rm gen}+\sum_{A}S_{A}=S_{\rm tot}$.

Generating the dynamical equations

Instances of the Community class can be initialized with any set of differential equations, which are specified as functions of the system state that return the time derivatives $dN_{i}/dt$ and $dR_{\alpha}/dt$. The package includes tools for constructing these functions automatically based on a dictionary of assumptions. These built-in dynamics are based on the recently introduced Microbial Consumer Resource Model (MicroCRM) Goldford2018 ; marsland2018available illustrated in Fig 2, which generalizes the classic consumer resource model of MacArthur and Levins MacArthur1970 to the microbial context by allowing organisms to release metabolic byproducts. Table 1 lists all the parameters of this family of models, along with the corresponding units.

In order to provide a general-purpose set of models that produce physically reasonable results, the MicroCRM assumes that all resource types are substitutable, and can all be converted to a common energy currency. This allows us to enforce energy conservation, preventing communities from bootstrapping themselves to large population sizes using metabolic secretions with no external resource supply. It also eliminates the need to specify in detail how each resource type interacts with all the others within the consumer metabolism. If such interactions are important for capturing a given experimental phenomenon, the built-in dynamics cannot be used, and custom functions must be written for $dN_{i}/dt$ and $dR_{\alpha}/dt$. The tutorial notebook included with the package contains an example of this kind, using Liebig’s Law of the Minimum to model phytoplankton dynamics.

Energy fluxes and growth rates

We begin by defining an energy flux into a cell $J^{\mathrm{in}}$, an energy flux that is used for growth $J^{\mathrm{growth}}$, and an outgoing energy flux due to byproduct secretion $J^{\mathrm{out}}$. Energy conservation requires

[TABLE]

for any reasonable metabolic model. It is useful to denote the input and output energy fluxes that are consumed/secreted in metabolite $\beta$ by $J_{\beta}^{\mathrm{in}}$ and $J_{\beta}^{\mathrm{out}}$ respectively. We can define corresponding mass fluxes by

[TABLE]

and

[TABLE]

where the conversion factor $w_{\beta}$ measures the energy density of metabolite $\beta$. In general, all these fluxes depend on the consumer species under consideration, and will carry an extra Roman index $i$ indicating the species.

We assume that a fixed quantity $m_{i}$ of power per cell is required for maintenance of species $i$, and that the per-capita growth rate is proportional to the remaining energy flux $(J^{\rm growth}-m_{i})$, with proportionality constant $g_{i}$. Under these assumptions, the time-evolution of the population size $N_{i}$ of species $i$ can be modeled using the equation

[TABLE]

We can model the resource dynamics by functions of the form

[TABLE]

where the function $h_{\alpha}$ describes the resource dynamics in the absence of consumers. The Community Simulator has two kinds of default resource dynamics: externally supplied and self-renewing. For externally supplied resources, we take a linearized form of the dynamics:

[TABLE]

while for self-renewing we take a logistic form

[TABLE]

Finally, the intrinsic resource dynamics can also be turned off, with $h_{\alpha}^{\rm off}=0$, to simulate resource depletion in a closed community with no resupply.

Input fluxes and output partitioning

We now specify the form of the input fluxes $\nu_{\beta}^{\mathrm{in}}$, and of the relationships among input, output and growth that define the metabolism. We start by assuming that all resource utilization pathways are independent, resulting in input fluxes of the form

[TABLE]

where $\sigma$ is a single-valued function encoding the relationship between resource availability and uptake rates. The community simulator implements three kinds of response functions: Type-I, linear response functions where

[TABLE]

a Type-II saturating Monod function,

[TABLE]

and a Type-III Hill or sigmoid-like function

[TABLE]

where $n>1$.

To obtain the output fluxes, we define a leakage fraction $l$ such that

[TABLE]

We allow different resources to have different leakage fractions $l_{\alpha}$. A direct consequence of energy conservation (Eq (1)) is that

[TABLE]

Finally, we denote by $D_{\beta\alpha}$ the fraction of the output energy that is contained in metabolite $\beta$ when a cell consumes $\alpha$. Note that by definition $\sum_{\beta}D_{\beta\alpha}=1$. The total energy output in metabolite $\beta$ is thus

[TABLE]

This also yields

[TABLE]

We are now in position to write down the full dynamics in terms of these quantities:

[TABLE]

Notice that when $\sigma$ is Type-I (linear) and $l_{\alpha}=0$ for all $\alpha$ (no leakage or byproducts), this reduces to MacArthur’s original model MacArthur1970 .

Metabolic regulation

The package can also generate dynamics for active metabolic regulation, which allocates a higher fraction of import capacity to nutrients with higher available energy flux. This regulation is implemented through a series of weight functions for resource $\alpha$ that reflect how much of the utilizable energy in the environment is in resource $\alpha$

[TABLE]

with $n_{\rm reg}$ a Hill coefficient that tunes steepness. Another option is to regulate based on the fraction of biomass contained in resource $\alpha$,

[TABLE]

For the metabolically regulated model, we define the input fluxes by

[TABLE]

Then, we can follow the exact same procedure as above. This yields the equations

[TABLE]

These equations are generated by the functions MakeConsumerDynamics and MakeResourceDyanamics, based on the user’s specification of the resource replenishment mode $h$, the response function $\sigma$, and the regulation mode $u$.

Sampling the parameters

The MicroCRM contains a large number of parameters: the $S_{\rm tot}$-dimensional vectors $g_{i}$ and $m_{i}$, the $M$-dimensional vectors $R_{\alpha}^{0},l_{\alpha},w_{\alpha}$ and $\tau_{\alpha}$ or $r_{\alpha}$, the $S_{\rm tot}\times M$ consumer preference matrix $c_{i\alpha}$ and the $M\times M$ metabolic matrix $D_{\alpha\beta}$. Some modeling choices require a small number of additional parameters: the maximal uptake rate $\sigma_{\rm max}$ for Type-II and Type-III growth, and the exponents $n$ for Type-III growth and $n_{\rm reg}$ for metabolic regulation. A dictionary containing all these parameters must be supplied to the Community instance upon initialization. A list of dictionaries may be supplied instead, to allow different wells to have different parameters.

The package contains a function MakeMatrices for generating the two matrices, which contain most of the ecological structure, based on a dictionary of modeling assumptions summarized in Table 2. The output of this function is illustrated in Fig 3 and described in detail below.

Consumer preferences $c_{i\alpha}$

We choose consumer preferences $c_{i\alpha}$ as follows. As stated earlier, we assume that each specialist family has a preference for one resource class $A$ (where $A=1\ldots F$) with $0\leq F\leq T$, and we denote the consumer coefficients for this family by $c_{i\alpha}^{A}$. We also consider generalists that have no preferences, with consumer coefficients $c_{i\alpha}^{\mathrm{gen}}$. The $c_{i\alpha}^{A}$ can be drawn from one of three probability distributions : (i) a Normal/Gaussian distribution, (ii) a Gamma distribution (which ensure positivity of the coefficients), and (iii) a Bernoulli distribution with binary preference levels. Fig 3 shows examples of all three models.

The Gaussian model is parameterized in terms of the mean $\mu_{c}=\left\langle\sum_{\alpha}c_{i\alpha}\right\rangle$ and variance $\sigma_{c}^{2}={\rm var}\left(\sum_{\alpha}c_{i\alpha}\right)$ of the total consumption capacity, and a parameter $q$ that controls how specialized each family is for its preferred resource class. In the generalist family, the mean and variance of $c_{i\alpha}$ are the same for all resources, and are given by

[TABLE]

where $\delta c_{i\alpha}^{\rm gen}=c_{i\alpha}^{\rm gen}-\langle c_{i\alpha}^{\mathrm{gen}}\rangle$ is the deviation from the mean value. The specialist families sample from a distribution with a larger mean for resources in their preferred class:

[TABLE]

where $M_{A}$ is the number of resources in class $A$ and $\mathbf{A}$ is the set of resource indices in class $A$. The variances are likewise larger for the preferred class:

[TABLE]

This makes it possible to construct pure specialist families with no off-target consumption by setting $q=1$. Note that this is different from the original version of this model in marsland2018available , where all the variances were chosen to be identical.

We also consider the case where consumer preferences are drawn from Gamma distributions, which guarantee that all coefficients are positive. Since the Gamma distribution only has two parameters, it is fully determined once the mean and variance are specified. We parameterize the mean and variance for this model in the same way as for the Gaussian model.

In the binary model, there are only two possible values for each $c_{i\alpha}$: a low level $\frac{c_{0}}{M}$ and a high level $\frac{c_{0}}{M}+c_{1}$. The elements of $c_{i\alpha}^{A}$ are given by

[TABLE]

where $X_{i\alpha}$ is a binary random variable that equals 1 with probability

[TABLE]

for the specialist families, and

[TABLE]

for the generalists. Note that the variance in each family is $\langle(\delta c_{i\alpha}^{A})^{2}\rangle=c_{1}^{2}p_{i\alpha}^{A}(1-p_{i\alpha}^{A})\sim c_{1}^{2}p_{i\alpha}^{A}$ for large $M$, which depends on $q$ in the same way as the variances in the Gaussian case.

Metabolic matrix $D_{\alpha\beta}$

We choose the metabolic matrix $D_{\alpha\beta}$ according to a three-tiered secretion model illustrated in Fig 3. The first tier is a preferred class of ‘waste’ products, such as carboyxlic acids for fermentative and respiro-fermentative bacteria, with $M_{w}$ members. The second tier contains byproducts of the same class as the input resource. For example, this could be attributed to the partial oxidation of sugars into sugar alcohols, or the antiporter behavior of various amino acid transporters. The third tier includes everything else. We encode this structure in $D_{\alpha\beta}$ by sampling each column $\beta$ of the matrix from a Dirichlet distribution with concentration parameters $d_{\alpha\beta}$ that depend on the byproduct tier, so that on average a fraction $f_{w}$ of the secreted flux goes to the first tier, while a fraction $f_{s}$ goes to the second tier, and the rest goes to the third. The Dirichlet distribution has the property that each sampled vector sums to 1, making it a natural way of randomly allocating a fixed total quantity (such as the total secretion flux from a given input). To write the expressions for these parameters explicitly, we let $A(\alpha)$ represent the class containing resource $\alpha$, and let $w$ represent the ‘waste’ class. We also introduce a parameter $s$ that controls the sparsity of the reaction network, ranging from a dense network with all-to-all connection when $s\to 0$, to maximal sparsity with each input resource having just one randomly chosen output resource as $s\to 1$. With this notation, we have

[TABLE]

The final two lines handle the case when the ‘waste’ type is being consumed. For these columns, the first and second tiers are identical. This led to an ambiguity in the expression presented in the Supporting Information of marsland2018available , which we have now clarified by treating this case separately. Note that in the third line, it is implicit that $A(\alpha)\neq w$, since $A(\alpha)=w$ is covered in the first line.

Propagation in time

Once an instance of the Community class is initialized, its state can be propagated forward in time using the Propagate method, as illustrated in Fig 1. Since the dynamical equations and parameters were supplied at initialization, the only required argument for this method is the time $T$. When the method is invoked, the state and parameters for each well are sent to different CPU’s (as many as are available) using the Pool.map function from the multiprocessing module in the Python standard library. Then the dynamical equations are integrated using the odeint function from SciPy, which calls the LSODA solver from the FORTRAN library ODEPACK scipy ; hindmarsh1983odepack .

If the initial population size of a species equals zero, but its per-capita growth rate is positive under current environmental conditions, the limited precision of any numerical solver will enable that species to spontaneously invade the community. To prevent this from happening, the Propagate method comes with an option compress_species (set to True by default), which reduces the dimensionality of the state and parameters before invoking the solver by removing all references to extinct species. Compression requires deciding whether each dimension of each parameter array corresponds to species or resources. This information is built in to the package for the MicroCRM, so c, D,w, g, m, l, R0, tau and r are all automatically handled correctly according to their definitions in that context. If custom dynamics are used, a dictionary of parameter dimensions must be supplied to the optional argument dimensions of Community when the instance is initialized, listing which parameters are of length S, length M, or of shape SxM, SxS, or MxM, respectively.

Passaging to fresh plates

The second method that can be invoked on a Community instance is Passage, which simulates pipetting of cultures to fresh wells in a typical 96-well plate experiment. The only required argument for Passage is a two-dimensional array f, whose elements $f_{\mu\nu}$ specify the fraction of the contents of well $\nu$ from the old plate that should be transfered to well $\mu$ on the new plate. The method also contains an option refresh_resource, set to True by default, that supplies the new plate with the same initial resource concentrations as the original one. This is the most direct way of simulating actual 96-well plate experiments, where the resources are resupplied in discrete intervals at each passaging step.

This method facilitates simulation of various kinds of mixing or coalescence experiments, as well as metacommunity dynamics of weakly coupled local communities. Fig 4 illustrates how this feature can capture coarse-grained spatial structure, following an experimental protocol developed for mimicking range expansions in 96-well plates datta2013range .

In addition, passaging helps to stabilize long simulations, and simulate demographic noise, by setting all cell counts to integer values. The species abundances $N_{i}$ are converted to absolute cell counts through a conversion factor scale, which is by default set to $10^{6}$. Then the new integer cell counts are generated by multinomial sampling based on the average number of cells $\sum_{\nu}f_{\mu\nu}(N_{i})_{\nu}$ of species $i$ transferred to well $\mu$. One source of numerical instability in ecological models is the exponentially small values reached by population sizes headed for extinction. The multinomial sampling ensures that these population sizes are fixed at zero when they become significantly smaller than 1 in absolute units. Thus a simple way of avoiding instability in a continuously resupplied chemostat model is to periodically call Passage with f set to the identity matrix and refresh_resource set to False.

Since the most common experiments involve many iterations of identical Passage and Propagate steps, the package also includes a method RunExperiment(f,T,np), which applies a given transfer matrix f and propagation time T for np iterations, saving a snapshot of the plate after each propagation step. If passaging is being used simply to stabilize the integration as discussed above, and not to simulate a batch culture setting, then the value of $T$ does not affect the results (as long as it is short enough to successfully eliminate instabilities). In this case, this variable mainly serves to control the time-resolution of the resulting timeseries.

Finding equilibrium points with convex optimization and expectation maximization

In many experimental contexts, one is interested in the stable community structure reached after a long period of constant environmental conditions. Numerical integration of the dynamical equations is an inefficient way to identify these equilibrium points. When the species diversity and number of resource types are high, the equations typically contain complex transient dynamics spanning a large range of time scales. This transient behavior is often irrelevant to the identification of the final equilibrium point, and wastes significant computation time. For a typical implementation of the MicroCRM with Type-I response and a $1,280\times 1,280$ binary consumer matrix, integrating to the steady state takes about 37 hours on standard hardware. The computation time appears to scale asymptotically as $M^{4}$ when the number of species $S$ and the number of resources $M$ are changed simultaneously, as shown in Fig. 6.

To address this problem, we have developed an algorithm for identifying equilibrium points directly, without integration through the transient. This algorithm is implemented in Community Simulator as the method SteadyState, and is illustrated in Fig. 5. Under the same test conditions, this algorithm converges between one and two orders of magnitude faster than numerical integration, as shown in Fig. 6, facilitating rapid hypothesis evaluation and iteration in large ecosystems. Note that while all other features of the Community Simulator depend only on packages included with a standard Anaconda installation (https://www.anaconda.com/), SteadyState additionally requires prior installation of a convex optimization package called CVXPY (https://www.cvxpy.org/).

The algorithm exploits a recently discovered duality between consumer resource models and constrained optimization over resource space mehta2018constrained ; Marsland2019a , which generalizes a minimization principle originally identified by MacArthur in the context of his original Consumer Resource Model MacArthur1970 . This duality applies to a wide class of consumer-resource type models, requiring only that the environmentally mediated interactions between pairs of consumer species are symmetric Marsland2019a .

For models in this class, it was shown that the vector of resource abundances $\mathbf{R}^{*}$ in every stable equilibrium state locally minimizes a measure $d(\mathbf{R}^{0},\mathbf{R})$ of the dissimilarity between the current resource abundances $\mathbf{R}$ and the supply point $\mathbf{R}^{0}$ (defined in general by $h_{\alpha}(R_{\alpha}^{0})=0$), subject to the constraint that all consumer growth rates are zero or negative ($dN_{i}/dt\leq 0$ for all $i$). Instances of the MicroCRM with Type I resource consumption, no metabolic regulation and $l_{\alpha}=0$ fall into this class. For externally supplied resources, $d$ turns out to be a weighted Kullback-Leibler (KL) divergence:

[TABLE]

while for self-renewing resources, it is a weighted Euclidean distance Marsland2019a :

[TABLE]

The equilibrium consumer populations $N_{i}^{*}$ are the Lagrange multipliers that enforce the constraints. For models with Type-I response, the non-invadable region is convex, allowing for efficient solution of the optimization problem using the Python package CVXPY cvxpy ; cvxpy_rewriting .

The duality does not strictly apply to other variants of the MicroCRM. In particular, byproduct secretion breaks the symmetry of the effective interactions between consumers whenever $l_{\alpha}>0$ for some resource $\alpha$ . The duality can be recovered, however, if the equilibrium point $\mathbf{R}^{0}$ of the intrinsic resource dynamics is changed to a new value $\tilde{\mathbf{R}}^{0}(\mathbf{R}^{*})$, which accounts for the extra resources produced by the consumer species when the system is at its equilibrium state $\mathbf{R}^{*}$ Marsland2019a . This accounting can be done in a variety of ways that all successfully recover the duality, with varying degrees of computational efficiency. The currently implemented form is

[TABLE]

where

[TABLE]

and $Q^{-1}_{\alpha\beta}$ are the elements of the matrix inverse of $Q_{\alpha\beta}$, satisfying $\sum_{\beta}Q^{-1}_{\alpha\beta}Q_{\beta\gamma}=\delta_{\alpha\gamma}$. With these definitions, one can show that the MicroCRM with Type I consumption and no regulation minimizes the objective function

[TABLE]

where

[TABLE]

To find $\mathbf{R}^{*}$, one must now self-consistently solve the following equation:

[TABLE]

The structure of this problem is mathematically equivalent to a standard task in machine learning, where one attempts to infer model parameters from partial data mehta2018high . These parameters $\theta$ specify a multivariate probability distribution $p(\mathbf{y}|\theta)$ for a set of measurements $\mathbf{y}$. A standard way of estimating the parameters is to compute the values $\hat{\theta}$ that maximize the likelihood of the data: $\hat{\theta}=\underset{\theta}{\rm argmax}\,p(\mathbf{y}|\theta)$. But if one actually has access to only a subset $\mathbf{x}$ of the measurement results, then the values $\mathbf{z}$ of the remaining quantities must also be estimated in order to perform this optimization. Ideally, one would use the statistical model with the optimal parameters $\hat{\theta}$ for this task. In the simplest case, where the value of $\mathbf{z}$ can be inferred with certainty given $\theta$ and $\mathbf{x}$, this results in the following self-consistency equation:

[TABLE]

This is identical in form to Eq. 36, where the parameters $\theta$ become the resource concentrations $\mathbf{R}$ and the estimated latent variables $\mathbf{z}$ become the effective unperturbed state $\tilde{\mathbf{R}}^{0}$. The observed data $\mathbf{x}$ become the model parameters, which are implicitly used in the calculation of $d$ and $\tilde{\mathbf{R}}$. We can think of the environmental perturbation $d$ as a statistical potential or “free energy” $-\ln p$, which is minimized when $p$ is maximized.

Eq. 37 can be solved by a standard iterative approach called Expectation Maximization mehta2018high . At each iteration $t$, the latent variable $\mathbf{z}_{t}$ is computed from the previous estimate $\hat{\theta}_{t-1}$ of $\hat{\theta}$, and then the new parameter estimate $\hat{\theta}_{t}$ is found by maximizing $p(\mathbf{x},\mathbf{z}_{t}|\theta)$. Fig. 5(c) contains pseudocode for this algorithm as applied to our ecological problem, which was also reported previously in Marsland2019a .

This algorithm fails to converge at low resource supply levels, because both arguments must be positive when $d$ is a weighted KL divergence, but $\tilde{R}^{0}_{t}$ can temporarily become negative under these conditions. To solve this issue, we replaced update step for $\tilde{R}^{0}_{t}$ by

  $\tilde{\mathbf{R}}^{0}_{t}\leftarrow\alpha\tilde{\mathbf{R}}^{0}(\mathbf{R}^{*}_{t})+(1-\alpha)\tilde{\mathbf{R}}^{0}_{t-1}$

where $\alpha$ is a constant rate, equivalent to the “learning rate” in machine learning mehta2018high .

Default values of the tolerance $\delta$ and learning rate $\alpha$ are set to $10^{-7}$ and 0.5, respectively, which give robust convergence for typical simulation scenarios. They can be adjusted as optional arguments of the SteadyState method.

This algorithm can be applied to any consumer-resource type model, including models beyond the MicroCRM framework, with non-substitutable resources Marsland2019a . But the enhanced efficiency of the new approach requires that the optimization problem be convex. In the Community Simulator package, the algorithm is only implemented for Type-I response with no metabolic regulation, where convexity is guaranteed. For more complex models, the differential equations must be numerically integrated using the Propagate method discussed above.

In the tutorial notebook included with the package, we show that the MicroCRM can be bistable if the externally supplied resources are insufficient to directly support growth of any consumer species. In this scenario the state with all consumers extinct is a stable equilibrium of the dynamics, and another stable equilibrium with persisting consumers is also possible that relies on the metabolic byproducts. In this scenario SteadyState method can find either of the two equilibria, depending on the initial estimate of $\tilde{\mathbf{R}}^{0}$, which can be set by an optional argument. If this initial condition is sufficiently close to the actual equilibrium state $\mathbf{R}^{0}$ of the intrinsic resource dynamics, the method ends in the state with the consumers extinct, where $\tilde{\mathbf{R}}^{0}=\mathbf{R}^{0}=\mathbf{R}^{*}$. But if the initial condition is not deliberately tuned to be close to $\mathbf{R}^{0}$, we find that the method typically finds the other state where some consumers survive.

Discussion

One interesting future direction to explore is integrating the Community Simulator with methods for directly analyzing Microbiome sequencing data. For example, there has been a renewed interest in statistical techniques such as Approximate Bayesian Computation (ABC) for understanding ecology and evolution csillery2010approximate . In ABC, the need to exactly calculate complicated likelihood functions – often a prerequisite for many statistical techniques – is replaced with the calculation of summary statistics and numerical simulations. For this reason, the Community Simulator Python package is ideally suited to form the backbone of new inference techniques for trying to related ecological processes to observed abundance patterns in microbial ecosystems.

Conclusion

We hope that the Community Simulator will become a valuable resource for the microbial ecology community. It has already played an important role in our own work. The package initially facilitated the systematic evaluation of the robustness of results to different modeling assumptions in a study of the effects of total energy influx on community structure, diversity and function Marsland2019a . More recently, the convex optimization approach has made it possible to perform more than 100,000 independent simulations in a reinterpretation and extension of Robert May’s classic work on diversity and stability Cui2019 ; May1972 . We have also employed the package to reproduce large-scale patterns in microbial biodiversity from the Human Microbiome Project, Earth Microbiome Project, and similar surveys marsland2019minimal . Finally, the random matrix approach implemented in this package is amenable to analytic calculation in the limit of large numbers of species and resources, using cavity methods from the physics of disordered systems Advani2018 ; cui2019effect . It is our belief that the Community Simulator will facilitate the further development of these mathematical techniques through efficient testing of new conjectures.

Availability and Requirements

• •

Project name: Community Simulator

• •

Project home page: https://github.com/Emergent-Behaviors-in-Biology/community-simulator

• •

Operating system(s): Linux or Mac preferred. Parallelization scheme is currently incompatible with Windows, and must be deactivated (set parallel=False when initializing a plate) for the code to run.

• •

Programming language: Python 3

• •

Other requirements: Numpy 1.15+, Pandas 0.23.0+, Matplotlib 2.2.3+, SciPy 1.1.0+. SteadyState method additionally requires CVXPY 1.0+.

• •

License: MIT

• •

Any restrictions to use by non-academics: None

Acknowledgements

This work was supported by NIH NIGMS grant 1R35GM119461 and Simons Investigator in the Mathematical Modeling of Living Systems (MMLS) award to PM. We are grateful to Kirill Korolev, Alvaro Sanchez, and Daniel Segrè for many useful conversations, and to Matti Gralka for testing cross-platform compatibility of the package. The performance evaluation reported in Fig 6 was performed on the Shared Computing Cluster which is administered by Boston University Research Computing Services.