Overview

MicroPop is an R package which simulates/predicts the growth of interacting populations of microbiota.

It is described in the paper:

Kettle H, G Holtrop, P Louis, HJ Flint. 2018. microPop: Modelling microbial populations and communities in R. Methods in Ecology and Evolution, 9(2), p399-409. doi: 10.1111/2041-210X.12873 https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12873

The user specifies the system via a number of input files (csv files that become dataframes) and the function microPopModel() will construct and solve the necessary equations (ordinary differential equations) and provide an output containing the solution (e.g. the concentrations of microbes, substrates and metabolites at the required time points) as well as all the settings/parameters involved in the simulation, and plots of the microbes and resources over time.

Getting started

Add the library to your current session

library(microPop)

#check package version 
packageVersion("microPop")
#> [1] '1.6'

Make sure you are using version 1.6 or above.

Defining your system

We will start with a very simple example.

We begin by making a csv file for our microbe ('MFG1.csv'; file provided in the InputFiles folder) and then loading it into R using the microPop function 'createDF()'. Note, it is better to use createDF() than the base R functions as no need to specify stringsAsFactors=FALSE etc.

M1=createDF(file='MFG1.csv')

M1 looks like this:

Where

• Rtype is the type of resource (S=substrate, P=product). Since products are often also substrates (e.g. in cross feeding) we use 'resources' to refer to both.

• keyResource is only important for specific cases (more later)

• pHcorners is only important if you want pH to limit growth (more later)

• numPathways is only important if the microbe has more than one metabolic pathway

• halfSat and maxGrowthRate are shown in plot below. The halfSat is the half saturation constant i.e. the amount of substrate needed to give half the maximum growth (assuming Monod Equation growth).

• The yield is g M1/g S1 i.e. the grams of substrate needed per gram of microbial growth.

• The units column is not used by microPop but is there for your reference.

Now we have specified the microbe we need to tell microPop about the environment it is in. We do this using 'system information (SI) files'. We can make these files in excel (or any editor) and save them as csv and read into R using createDF() again:

res.SI=createDF(file='resourceSI.csv')

res.SI looks like this:

Note, we have specified the starting value of each resource, the amount that is flowing into the system and the specific wash out rate (e.g. if the system has a turnover time of 2 days then this is 0.5 /d and each resource, R, (g/l) leaves at the rate 0.5*R g/l/d ). For convenience we also specify the molar mass of each resource here.

We now read in the microbial SI file:

mic.SI=createDF(file="microbeSI.csv")

mic.SI looks like:

Note, once these input files are defined as objects in our workspace we can edit them any way we like in R.

We now use the input arguments in the microPopModel() function to define the remaining information about the system i.e. which microbes are there (microbeNames) and how long we want to simulate (times) and we tell it which SI data frames to use. Note, here times is a sequence from 0 to 1 going up in 0.001 intervals. There are no units specified but since our input files all have time units of d then this is what microPop assumes are the units for times. You can use any time units you like but you must be consistent throughout all your input information. Similarly only use one mass unit. MicroPop is based on mass balance equations therefore you can not use moles only mass or concentrations for quantities.

out=microPopModel(microbeNames = 'M1',
              times=seq(0,1,0.001),
              resourceSysInfo=res.SI,
              microbeSysInfo=mic.SI,
              plotOptions=list(plotFig=FALSE)
)

par(mfrow=c(1,2))
plotMicrobes(out,cex.title=0.7,cex.legend=0.5,cex.ax=0.7)
plotResources(out,cex.title=0.7,cex.legend=0.5,cex.ax=0.7)

When microPop runs it will automatically produce two plots - one for the resources (i.e. substrates and products) and another for the microbes. This does not work well with R markdown (used to contruct this vignette) so we turn this off using plotOptions(plotFig=FALSE) and then call the two plotting functions ( plotMicrobes and plotResources) separately afterwards (note, we will not always show this plotting code in the vignette from now on).

Let's look at the output ('out'):

names(out)
#> [1] "solution" "parms"

It is a list with 2 elements. The first item, 'solution', is a matrix of the state variables (M1, S1 and P1) at the times specified in 'times'. Let's look at the first few lines:

head(out$solution)
#>       time       M1        S1         P1 pH
#> [1,] 0.000 1.000000 10.000000 0.00000000 NA
#> [2,] 0.001 1.004008  9.974314 0.01168303 NA
#> [3,] 0.002 1.008031  9.948588 0.02340120 NA
#> [4,] 0.003 1.012071  9.922820 0.03515467 NA
#> [5,] 0.004 1.016127  9.897010 0.04694356 NA
#> [6,] 0.005 1.020199  9.871158 0.05876805 NA

If you are familiar with R you can make your own plots from the matrix out$solution, e.g.

plot(range(out$solution[,'time']),c(0,10),type='n',xlab='Time (d)',ylab='concentration (g/l)')
lines(out$solution[,'time'],out$solution[,'M1'],col='black')
lines(out$solution[,'time'],out$solution[,'S1'],col='red')
lines(out$solution[,'time'],out$solution[,'P1'],col='blue')
legend('topright',c('M1','S1','P1'),col=c('black','red','blue'),lty=1)

The second item, 'parms' is another list containing all the information about the simulation. We can look at the names of its entries:

names(out$parms)
#>  [1] "keyRes"               "numPaths"             "allStrainNames"      
#>  [4] "resourceNames"        "microbeNames"         "numStrains"          
#>  [7] "Pmats"                "Smats"                "strainPHcorners"     
#> [10] "nonBoostFrac"         "pHLimit"              "pHVal"               
#> [13] "molarMass"            "bacCutOff"            "balanceTol"          
#> [16] "checkMassConv"        "times"                "networkAnalysis"     
#> [19] "myPars"               "pHFunc"               "pHLimFunc"           
#> [22] "extraGrowthLimFunc"   "growthLimFunc"        "combineGrowthLimFunc"
#> [25] "uptakeFunc"           "productionFunc"       "combinePathsFunc"    
#> [28] "massBalanceFunc"      "entryRateFunc"        "removalRateFunc"

You can look at any of these entries to check the model settings for the simulation, e.g.

out$parms$microbeNames
#> [1] "M1"
out$parms$resourceNames
#> [1] "S1" "P1"

Let's see what happens if we add another microbe to this system.

The new microbe (M2) has a high maximum growth rate but also a high half-saturation constant compared with M1.

M2=createDF(file="MFG2.csv")

M2 looks like this:

We also need to add this microbe to the microbe Sys Info file. Let's assume for simplicity that it has the same conditions as M1, then we can add the M1 column to the dataframe and name it M2:

micNew.SI = cbind(mic.SI,M2=mic.SI[,'M1']) #cbind is a function which binds columns

We now run a simulation with the two microbes and extend the time period to 3 days

out=microPopModel(microbeNames = c('M1','M2'),
              times=seq(0,3,0.001),
              resourceSysInfo=res.SI,
              microbeSysInfo=micNew.SI,
              plotOptions = list(plotFig=FALSE)
)

#SOLUTION 
#look at the final concentrations:

out$solution[nrow(out$solution),c('M1','M2')]
#>        M1        M2 
#> 0.3589012 0.1750382

#M1 and M2 have the same yield but M2 has a higher max growth rate and a higher halfSat. 
#As time goes on the substrate concentration is low so the max growth rate is not important.
#We change half-sat of M2 to 1e-4 (i.e. less than M1) and look what happens:
newM2=M2 #make a new copy of M2 to play with
newM2['halfSat','S1']=1e-4
micNew.SI = cbind(micNew.SI,newM2=micNew.SI[,'M2']) #add newM2 to the microbe SI data frame (assume same conditions as M2)

out=microPopModel(microbeNames = c('M1','newM2'),
                  times=seq(0,3,0.001),
                  resourceSysInfo=res.SI,
                  microbeSysInfo=micNew.SI,
                  plotOptions = list(plotFig=FALSE)
)
#> [1] "Set up completed, ODE solver called..."

out$solution[nrow(out$solution),c('M1','newM2')]
#>         M1      newM2 
#> 0.07041673 0.46358173
#Now M2 > M1 near the end of the simulation

#SOLUTION 
out=microPopModel(microbeNames = c('M1','M2'),
                  times=seq(0,10,0.001),
                  resourceSysInfo=res.SI,
                  microbeSysInfo=micNew.SI,
                  plotOptions = list(plotFig=FALSE)
)
#> [1] "Set up completed, ODE solver called..."

out$solution[nrow(out$solution),c('M1','M2')]
#>           M1           M2 
#> 0.2999732783 0.0001651169
#M2 dies out as it can't compete with M1 at low substrate concentrations.

#SOLUTION
#start by making M1 and M2 the same:
newM2[c('halfSat','yield','maxGrowthRate'),'S1'] = M1[c('halfSat','yield','maxGrowthRate'),'S1'] 

out=microPopModel(microbeNames = c('M1','newM2'),
                  times=seq(0,6,0.001),
                  resourceSysInfo=res.SI,
                  microbeSysInfo=micNew.SI,
                  plotOptions=list(plotFig=FALSE)
)
#> [1] "Set up completed, ODE solver called..."

#Look at steady state values
out$solution[nrow(out$solution),c('M1','newM2','S1')]
#>           M1        newM2           S1 
#> 0.1557893117 0.1557893117 0.0002384204


#Now investigate the effect of yield values:

#make M2 have a lower yield than M1
newM2['yield','S1'] = 0.5*as.numeric(M1['yield','S1'])

out=microPopModel(microbeNames = c('M1','newM2'),
                  times=seq(0,6,0.001),
                  resourceSysInfo=res.SI,
                  microbeSysInfo=micNew.SI,
                  plotOptions = list(plotFig=FALSE)
)
#> [1] "Set up completed, ODE solver called..."
out$solution[nrow(out$solution),c('M1','newM2','S1')]
#>           M1        newM2           S1 
#> 0.1046860252 0.1046860252 0.0002360937

Multiple Strains

When we define the input file for a microbe, this can be viewed as a single strain or as a functional group. If it is a functional group then we can add a number of strains to it. We do this by setting the numStrains argument in microPopModel() to more than 1. We can then specify the parameters, e.g. halfSat, yield, maxGrowthRate etc, for these strains to vary randomly around the average value for the group. This way we have a population of strains with the same metabolic paths but with slightly different growth traits.

out=microPopModel(microbeNames = c('M1','M2'),
              times=seq(0,3,0.001),
              resourceSysInfo=res.SI,
              microbeSysInfo=micNew.SI,
              numStrains=3,
              strainOptions=list(randomParams=c('halfSat','yield','maxGrowthRate'),
                                 seed=1,distribution='uniform',
                                 percentTraitRange=20,
                                 applyTraitTradeOffs=TRUE,
                                 tradeOffParams=c('halfSat','maxGrowthRate')),
              plotOptions = list(plotFig=FALSE)
)

par(mfrow=c(1,2))
plotMicrobes(out,cex.title=0.7,cex.legend=0.5,cex.ax=0.7,sumOverStrains = FALSE)
plotResources(out,cex.title=0.7,cex.legend=0.5,cex.ax=0.7)

Here there are 3 strains in each functional group (M1, M2) whose parameter values vary by +/- 20% from a uniform distribution around the mean value specified in the input data frames. The stochastically-generated values have been picked so there are trade-offs between the maximum growth rate and the half-saturation constant. Note that we have to specify sumOverStrains=FALSE to plot each strain in plotMicrobes and if we were using the plotting function run in _microPopModel we would use plotOptions = list(plotFig=TRUE,sumOverStrains=FALSE).

We can look at the final distribution of biomass:

barplot(out$solution[nrow(out$solution),out$parms$allStrainNames],col=c(rep('red',3),rep('cyan',3)))

Note the microPop nomenclature for naming the stochastic strains and these names are found in

out$parms$allStrainNames
#> [1] "M1.1" "M1.2" "M1.3" "M2.1" "M2.2" "M2.3"

It is also possible to have different numbers of strains in each group. In this case numStrains is a named vector,

out=microPopModel(microbeNames = c('M1','M2'),
              times=seq(0,3,0.001),
              resourceSysInfo=res.SI,
              microbeSysInfo=micNew.SI,
              numStrains=c(M1=2,M2=3),
              strainOptions=list(randomParams=c('halfSat','yield','maxGrowthRate'),
                                 seed=1,distribution='uniform',
                                 percentTraitRange=20,
                                 applyTraitTradeOffs=TRUE,
                                 tradeOffParams=c('halfSat','maxGrowthRate')),
              plotOptions = list(plotFig=FALSE)
)

#SOLUTION
out3=out #save the 3 day solution

out10=microPopModel(microbeNames = c('M1','M2'),
              times=seq(0,10,0.001),
              resourceSysInfo=res.SI,
              microbeSysInfo=micNew.SI,
              numStrains=3,
              strainOptions=list(randomParams=c('halfSat','yield','maxGrowthRate'),
                                 seed=1,distribution='uniform',
                                 percentTraitRange=20,
                                 applyTraitTradeOffs=TRUE,
                                 tradeOffParams=c('halfSat','maxGrowthRate')),
              plotOptions = list(plotFig=FALSE)
)
#> [1] "Set up completed, ODE solver called..."

barplot(out10$solution[nrow(out10$solution),out10$parms$allStrainNames],col=c(rep('red',3),rep('cyan',3)))

out20=microPopModel(microbeNames = c('M1','M2'),
              times=seq(0,20,0.001),
              resourceSysInfo=res.SI,
              microbeSysInfo=micNew.SI,
              numStrains=3,
              strainOptions=list(randomParams=c('halfSat','yield','maxGrowthRate'),
                                 seed=1,distribution='uniform',
                                 percentTraitRange=20,
                                 applyTraitTradeOffs=TRUE,
                                 tradeOffParams=c('halfSat','maxGrowthRate')),
              plotOptions = list(plotFig=FALSE)
)
#> [1] "Set up completed, ODE solver called..."

barplot(out20$solution[nrow(out20$solution),out20$parms$allStrainNames],col=c(rep('red',3),rep('cyan',3)))

#Plot all together
barplot(rbind(out3$solution[nrow(out3$solution),out3$parms$allStrainNames],
              out10$solution[nrow(out10$solution),out10$parms$allStrainNames],
              out20$solution[nrow(out20$solution),out20$parms$allStrainNames]),
        beside=TRUE,legend.text = c('3 days','10 days','20 days'))
#> Warning in rbind(out3$solution[nrow(out3$solution),
#> out3$parms$allStrainNames], : number of columns of result is not a multiple of
#> vector length (arg 1)

Rate Functions

MicroPop uses ordinary differential equations which define the rates of changes of the microbes, substrates and metabolic products (the "state variables") involved in the system. However, the general user does not need to concern themselves with these unless they want to change the default functions.

There are default functions for

• rate of entry of each state variable to the system -- entryRateFunc()

• rate of exit of each state variable -- removalRateFunc()

• the limit on the maximum growth rate e.g. Monod equation -- growthLimFunc()

• how to combine growth on multiple resources -- combineGrowthLimFunc()

• resource uptake rate due to microbial growth -- uptakeFunc()

• production rate of metabolites -- productionFunc()

• pH at each time point -- pHFunc()

• the pH limit on growth -- pHLimFunc()

• other limits on growth (e.g. inhibition) -- extraGrowthLimFunc()

• combining growth on multiple metabolic paths -- combinePathsFunc()

These are all contained in a list called rateFuncs which has a default version called rateFuncsDefault. You can see this if you type rateFuncsDefault at the prompt (it is quite long though!). To just see the elements you can type

  names(rateFuncsDefault)
#>  [1] "pHFunc"               "pHLimFunc"            "extraGrowthLimFunc"  
#>  [4] "growthLimFunc"        "combineGrowthLimFunc" "uptakeFunc"          
#>  [7] "productionFunc"       "combinePathsFunc"     "massBalanceFunc"     
#> [10] "entryRateFunc"        "removalRateFunc"

To look at one function, e.g. growthLimFunc use

  rateFuncsDefault$growthLimFunc
#> function (strainName, groupName, pathName, varName, resourceValues, 
#>     allSubType, strainHalfSat, stateVarValues, parms) 
#> {
#>     if (resourceValues[varName] <= 0) {
#>         v = 0
#>     }
#>     else {
#>         if (allSubType[varName] == "Sb" | allSubType[varName] == 
#>             "Se") {
#>             v = resourceValues[varName]/(resourceValues[varName] + 
#>                 strainHalfSat[varName])
#>         }
#>         if (allSubType[varName] == "S" | allSubType[varName] == 
#>             "Sm") {
#>             v = resourceValues[varName]/(strainHalfSat[varName] * 
#>                 (1 + sum(resourceValues[allSubType == "S" | allSubType == 
#>                   "Sm"]/strainHalfSat[allSubType == "S" | allSubType == 
#>                   "Sm"])))
#>         }
#>     }
#>     return(max(v, 0))
#> }
#> <bytecode: 0x77ec158>
#> <environment: namespace:microPop>

This looks complicated because microPop is very generic so there are different equations for dealing with different types of resource, e.g. substitutable (S) and essential (Se) resources. Also there may be boosting resource (Sb) e.g. Acetate, or for viruses a microbe may be a resource (Sm). However, should you wish to replace this with your own function (which must have the same form of input arguments and output) then this is done as follows:

  myRateFuncs = rateFuncsDefault #make your own copy of the rateFuncs list
  myRateFuncs$growthLimFunc = myGrowthFunction #assign your own function  

where myGrowthFunction is your own function (which you will need to define before this point). Then when you call microPop you use the new rateFuncs list , i.e.

out=microPopModel(microbeNames = 'M1',
              times=seq(0,10,0.001),
              resourceSysInfo=res.SI,
              microbeSysInfo=mic.SI,
              rateFuncs=myRateFuncs
)

pH

It is well established that pH has a significant effect on microbial growth rates. We now import a new microbe, M3, which has pH preferences, consumes 2 substrates and generates 2 products.

M3=createDF(file='MFG3.csv')

M3 looks like this:

The pH preference is defined by pHcorners which describe the following shape:

plot(as.numeric(M3['pHcorners',2:5]),c(0,1,1,0),type='l',xlab='pH',ylab='pH limitation',col='blue',lwd=2)

This multiplies the maxGrowthRate, thus Where the limitation is 1 there is no restriction on growth and where the limit is 0 there is no growth; inbetween it it linearly scaled.

We also need to read in the new resource SI file with the system information about these new resources (S2,P2).

res.SI2=createDF(file='resourceSI2.csv')

res.SI2:

mic.SI2=createDF(file='microbeSI2.csv')

mic.SI2

In order for the pH preferences to be relevant we need to specify the pH of the system. This can be done using the pHVal argument in microPopModel() (see code below) or if you want pH to change over time then you can use pHFunc in the rateFuncs list. To make sure the pH affects the growth we also need to set the microPopModel argument, pHLimit=TRUE. We now run simulations for different pHs and plot the results

out1=microPopModel(microbeNames = 'M3',
              times=seq(0,5,0.001),
              resourceSysInfo=res.SI2,
              microbeSysInfo=mic.SI2,
              pHLimit=TRUE,
              pHVal=4.1,
              plotOptions=list(plotFig=FALSE)
              )
    
out2=microPopModel(microbeNames = 'M3',
              times=seq(0,5,0.001),
              resourceSysInfo=res.SI2,
              microbeSysInfo=mic.SI2,
              pHLimit=TRUE,
              pHVal=6,
              plotOptions=list(plotFig=FALSE)
              )

plot(range(out1$solution[,'time']),c(0,max(c(out1$solution[,'M3'],out2$solution[,'M3']))),type='n',xlab='Time (d)',ylab='Concentration (g/l)', main='M3',lwd=2)
lines(out1$solution[,'time'],out1$solution[,'M3'],col='black',lwd=2)
lines(out2$solution[,'time'],out2$solution[,'M3'],col='red',lwd=2)
legend('right',col=1:2,lty=1,legend=c(4.1,6),title='pH',lwd=2)

#SOLUTION
myRateFuncs=rateFuncsDefault

myRateFuncs$pHFunc=function(time,parms){
  if (time<=2){
    pH=6
  }else{
    pH=4.1
  }
  return(pH)
}

out=microPopModel(microbeNames = 'M3',
              times=seq(0,5,0.001),
              resourceSysInfo=res.SI2,
              microbeSysInfo=mic.SI2,
              pHLimit=TRUE,
              rateFuncs=myRateFuncs
)
#> path1 
#>    "" 
#> [1] "MICROPOP NOTE TO USER: The key resource for M3 on path 1 is not used"
#> [1] "Set up completed, ODE solver called..."

Included Data and Demo files

The microPop package comes with many microbial data frames and system information files included. These are for modelling the human colonic microbiota, the rumen etc (they are described in Kettle et al, 2018 and used in the demofiles). You can get a list of the included data sets using:

data(package='microPop')

Or on the Environment tab in Rstudio, change the drop down menu from Global Environment to package:microPop and then you will see the data sets in the package. You can click on these to look at them. If you are not in Rstudio you can use View() to look at them or just type the name of the data frame at the R prompt.

There are also a number of demo files included in the package. These correspond to the examples in Kettle et al (2018). To find where they are on your computer, use

system.file('DemoFiles/',package='microPop')

You can then view them in any editor.

To run any of the examples in Kettle et al 2018, you can use the function runMicroPopExample(). E.g.

runMicroPopExample('human1') 

out=microPopModel(
    microbeNames=c('Bacteroides','NoButyStarchDeg','Acetogens'),
    times=seq(0,4,1/24),
    resourceSysInfo=resourceSysInfoHuman,
    microbeSysInfo=microbeSysInfoHuman,
    plotOptions = list(plotFig=FALSE)
)

simulation.times=seq(0,4,1/24)

myRateFuncs=rateFuncsDefault

myRateFuncs$pHFunc=function(time,parms){
    if (time<=max(simulation.times)/2){pH=5.5}else{pH=6.5}
    return(pH)
}

out=microPopModel(
    microbeNames=c('Bacteroides','NoButyStarchDeg','Acetogens'),
    times=simulation.times,
    resourceSysInfo=resourceSysInfoHuman,
    microbeSysInfo=microbeSysInfoHuman,
    rateFuncs=myRateFuncs,
    pHLimit=TRUE,
    plotOptions = list(plotFig=FALSE)
)

#SOLUTION

myBacteroides=Bacteroides

print(myBacteroides['pHcorners',])
#>           units Protein  NSP   RS Acetate Propionate Succinate H2 CO2 other
#> pHcorners    pH     5.6 6.35 7.85     8.6
print(myBacteroides['pHcorners',2:5])
#>           Protein  NSP   RS Acetate
#> pHcorners     5.6 6.35 7.85     8.6
myBacteroides['pHcorners',2:5]=c(5,6,7,8)

microbeHuman.SI = cbind(microbeSysInfoHuman,myBacteroides=microbeSysInfoHuman[,'Bacteroides']) #add Sys Info

out=microPopModel(
    microbeNames=c('myBacteroides','NoButyStarchDeg','Acetogens'),
    times=simulation.times,
    resourceSysInfo=resourceSysInfoHuman,
    microbeSysInfo=microbeHuman.SI,
    rateFuncs=myRateFuncs,
    pHLimit=TRUE
)
#> [1] "Set up completed, ODE solver called..."

Stoichiometries and essential resources

In the M3 dataframe there are two substrates S1, S2 and two products P1, P2. The molar masses of these are defined in res.SI2 and are S1=100 g/mol, S2=100 g/mol, P1=50 g/mol and P2=25 g/mol. The stoichiometry will only balance if either S1 or S2 is used - not both.

2 S1 or 2 S2 -> 3 P1 +2 P2

2 x 100g or 2 x 100g -> 3 x 50g + 2 x 25g

This is because these are substitutable resources.

The substrates are given Rtype S for substitutable, i.e. M3, can grow on either/both S1 and S2. This means that the stoichiom entry means 2 mol of S1 or 2 mol of S2 will give 3 mol of P1 and 2 mol of P2.

If both substrates were required for growth these would be given an Rtype of Se for essential substrate. E.g., methanogens require both CO2 and H2 to produce methane. Methanogens are an included data frame and can be viewed without loading. There is a version for the human colon (Methanogens) and for the rumen (Xh2). We will look at Xh2 as it is a simpler version (no formate path):

Notice here that the microbial biomass is also included as a product in the stoichiometry (and therefore ammonia is also included as a substrate). The prefix 'S' indicates a soluble element.

When we are dealing with essential resources the keyResource category is used by microPop. This is needed as it is the substrate used for the yield value i.e. if H2 (Sh2) is the keyResource then a yield of 0.0904 is 0.09g of methanogens per 1 g of H2.