System analysis: Competition and coexistence

Author

David Kneis (firstname.lastname @ tu-dresden.de)

Published

February 21, 2026

1 Introduction

According to classical theory, the key to co-existence is the presence of contrasting ecological niches. But even in seemingly simply-structured environments like the pelagic zone, we find substantial species diversity.

Does that diversity actually reflect stable co-existence between species? Or is it rather a case of mere co-presence of unequal competitors? In fact, this question is very hard to answer with field studies. First of all, investigation time is usually limited (no time to study long-term competition, unless its about microbes). Second, field observations are always contaminated with noise, reflecting various random events.

Therefore, models are a quite useful means to study the possible mechanisms and necessary conditions for co-existence. In particular, simulations can easily be run over very long periods of time such that even very minor selective advantages lead to ultimate outcompetition.

In spite of the apparent simplicity of models, the interpretation of results may be challenging! In this lecture, we focus on a small selection of simple cases and scenarios. This selection is not exhaustive but a reasonable starting point for understanding more complex situations.

To make the topic digestible at all, we build on technology for ODE-based modeling introduced earlier, namely the rodeo package.

2 General study system

We consider a continuous flow stirred tank reactor (CFSTR) again, so as to avoid all unnecessary complexity. In contrast to earlier lectures, we now consider two species sharing space and the single unique resource pumped into the reactor.

In a series of scenarios, we play with the particular capabilities of the species and the external forcings to see if and when stable coexistence can be achieved.

3 Technical basis of subsequent analyses

3.1 Concept

It is assumed that the model specification is in a workbook in .xlsx format with the two worksheets declarations and equations such that it can be imported by the buildFromWorkbook function of the rodeo package.

Two make the study of scenarios more convenient, I wrote the two wrapper functions below. They can be used to instantly run a dynamic simulation (dyn()) or steady-state computation (std()), including visualization, with minimum inputs. Nevertheless, the functions are flexible enough to allow parameters or initial values to be adjusted using their update argument.

library("rodeo")

# runs a dynamic simulation
dyn <- function(
  file,                       # workbook defining the ODE model
  update=NULL,                # named vector to alter parameters, initial v.
  times=seq(0, 240, 1),       # requested prediction times
  graphics=TRUE,              # shall results be plotted?
  maxcols=3,                  # max. columns in plot layout
  trans=function(x){x^0.25}   # function transforming y-axis in plots
) {
  m <- buildFromWorkbook(file)        # load model, including default values
  if (!is.null(update)) {             # possibly overwrite defaults
    if (is.null(names(update))) 
      stop("cannot update items without names")
    if (!all(names(update) %in% c(m$namesVars(), m$namesPars()))) 
      stop("cannot update value(s) of unknown items")
    change <- function(m, x, pars=TRUE) {
      what <- if (pars) "Pars" else "Vars"
      x <- x[names(x) %in% do.call(m[[paste0("names",what)]], list())]
      if (length(x) > 0) {
        z <- do.call(m[[paste0("get",what)]], list())
        z[names(x)] <- x
        do.call(m[[paste0("set",what)]], list(z))
      }
      m
    }
    m <- change(m, update[m$namesVars()], pars=FALSE)
    m <- change(m, update[m$namesPars()], pars=TRUE)
  }
  x <- m$dynamics(times)              # run dynamic simulation
  if (graphics) {
    nc <- min(m$lenVars(), maxcols)
    nr <- ceiling(m$lenVars() / nc)
    par(mfrow=c(nr,nc))              # plot with custom layout
    par(cex=1)
    for (v in m$namesVars()) {
      plot(x[,"time"], trans(x[,v]), type="l", bty="L",
        yaxt="n", xlab="time", ylab=v)
      yt <- c(0,10^(-6:9))
      axis(2, at=trans(yt), labels=yt)
    }
    par(mfrow=c(1,1))
    invisible(NULL)
  } else {
    x[nrow(x), m$namesVars()]         # return final state
  }
}

# computes a steady-state solution
std <- function(
  file,           # workbook defining the ODE model
  update=NULL,    # named vector to alter parameters, initial values
  tini=5000,      # assumed time to steady state tried initially
  tmul=4,         # factor applied to tmax in every new iteration
  abstol=1e6,     # stopping criterion, max. abs. change in any state v.
  maxiter=10      # max. number of iterations (emergency exit)
) {
  t <- tini
  i <- 0
  old <- dyn(file, update, times=c(0, t/2), graphics=FALSE)
  while (i < maxiter) {
    new <- dyn(file, update, times=c(0, t), graphics=FALSE)                   
    if (max(abs(old - new)) < abstol) {
      break
    } else {
      t <- t * tmul
      old <- new
      i <- i + 1
    }
  }
  if (i == maxiter) stop("no steady state found in 'maxiter' iterations")
  new
}

To make use of these functions in your own R script, just copy the code and paste it into an empty text file. Let this file be called utility_functions.r, for example. Then, you would make the functions available by including the following line at the top of your script like so:

source("utility_functions.r")

3.2 Demonstration of the `dyn` and `std` function

For the purpose of demonstration, we consider a single-species case again: a hypothetical species A is growing on a resource R. Growth of the species and flushing of the reactor are the only processes considered (Figure 1).

Table 1: State variables and parameters of the single species bioreactor model.

type	name	unit	description	default
variable	R	mg / mL	resource concentration	0.0e+00
variable	A	cells / mL	abundance of organisms	1.0e+05
parameter	mu	1 / h	max. growth rate	2.0e+00
parameter	h	mg / mL	half saturation const.	1.8e-01
parameter	y	cells / mg	yield of A growing on R	1.0e+06
parameter	q	mL / h	flow rate	1.0e+01
parameter	v	mL	reactor volume	1.0e+02
parameter	Rin	mg / mL	input resource conc.	1.5e+01
parameter	Ain	cells / mL	abundance of A in input	0.0e+00

Let the model be available as a workbook singleSpecies.xlsx with the contents of Table 1 in sheet declarations and the contents of Table 2 on sheet equations.

Table 2: Simultaneous differential equations of the single species bioreactor model.

name	unit	description	rate	A	R
flow	1 / h	im-/export	q / v	Ain-A	Rin - R
growth	cells / mL / h	growth	mu * R / (R + h) * A	1	-1 / y

Figure 1: Symbolic stoichiometry matrix of the single species bioreactor model.

Then, we can trigger a dynamic simulation, including visualization (Figure 2) with this single line of code:

dyn("models/singleSpecies.xlsx")

Figure 2: Prediction of the single-species model using default parameters and initial values.

To rerun the simulation with certain parameters or initial values being changed, we use the same call but specify the update argument. Here, we modify the flow rate and the growth rate at the same time.

dyn("models/singleSpecies.xlsx", update=c(q=100, mu=1.5))

Figure 3: Prediction of the single-species model with changed parameter values.

Finally, we can directly ask for the steady-state conditions (optionally making use of the update argument again).

std("models/singleSpecies.xlsx")

           R            A 
9.473684e-03 1.499053e+07

4 Scenarios with two competitors

4.1 Base scenario: Straight outcompetition

This scenario differs from the single species bioreactor (species A only) by the mere presence of a competitor species B. With regard to the considered processes, not much has changed: Suspended cells grow and they get flushed by the continuous flow (Figure 4).

Figure 4: Symbolic stoichiometry matrix of the basic two species bioreactor model.

Even the default parameters and initial values for A and B were chosen to be identical (Table 3, Table 4).

Table 3: State variables and parameters of the basic two species bioreactor model.

type	name	unit	description	default
variable	R	mg / mL	resource concentration	0.0e+00
variable	A	cells / mL	abundance, organism A	1.0e+05
variable	B	cells / mL	abundance, organism B	1.0e+05
parameter	mu_A	1 / h	max. growth rate of A	2.0e+00
parameter	mu_B	1 / h	max. growth rate of B	2.0e+00
parameter	h_A	mg / mL	half saturation const. of A	1.8e-01
parameter	h_B	mg / mL	half saturation const. of B	1.8e-01
parameter	y_A	cells / mg	yield of A on R	1.0e+06
parameter	y_B	cells / mg	yield of B on R	1.0e+06
parameter	q	mL / h	flow rate	1.0e+01
parameter	v	mL	reactor volume	1.0e+02
parameter	Ain	cells / mL	abundance of A in input	0.0e+00
parameter	Bin	cells / mL	abundance of B in input	0.0e+00
parameter	Rin	mg / mL	input resource conc.	1.5e+01

Table 4: Simultaneous differential equations of the basic two species bioreactor model.

name	unit	description	rate	A	B	R
flow	1 / h	im-/export	q / v	Ain-A	Bin-B	Rin - R
growth_A	cells / mL / h	growth	mu_A * R / (R + h_A) * A	1		-1 / y_A
growth_B	cells / mL / h	growth	mu_B * R / (R + h_B) * B		1	-1 / y_B

Running a dynamic simulation suggest that both competitors survive and the abundances of A and B in steady state are equal (Figure 5). This is not very surprising because, with identical default parameters, A and B are essentially indistinguishable as if they were the same species.

dyn("models/twoSpecies_basic.xlsx")

Figure 5: Prediction of the two-species model with defaults.

Now, let species A have a fitness advantage over B. We could create this situation, for example, by setting mu_A greater than mu_B. Alternatively, we can adjust the half saturation constant of the Monod model (h_B) such that B requires a higher resource concentration than A to reach an identical effective growth rate. Apparently, the superiority of A with regard to resource dependency drives B straight into extinction (Figure 6).

dyn("models/twoSpecies_basic.xlsx", update=c(h_B=0.36))

Figure 6: Prediction of the two-species model after compromising the competitive ability of species B.

This is easily verified by computing the steady state (note the very low value for species B).

std("models/twoSpecies_basic.xlsx", update=c(h_B=0.36))

           R            A            B 
9.473684e-03 1.499053e+07 1.156247e-12

Can we possibly compensate for the disadvantage of B by a higher inoculum? We can test this by modifying the initial proportions of A and B accordingly and rerun the model (Figure 6).

dyn("models/twoSpecies_basic.xlsx", update=c(h_B=0.36, A=1e3, B=1e8))

Figure 7: Prediction of the two-species model after compromising the competitive ability of species B and reducing the inoculumn of species A.

Apparently, the higher initial value does not save the survival of species B in the long run, as confirmed by the calculation of the steady state. The results do not change, no matter how extreme the initial proportion is chosen.

std("models/twoSpecies_basic.xlsx", update=c(h_B=0.36, A=1e3, B=1e8))

           R            A            B 
9.473684e-03 1.499053e+07 2.992268e-11

4.2 Forced co-presence by invasion

The simplest option to keep the inferior competitor in the system is to allow for a certain rate of invasion. In the model, we can simply mimic invasion by setting the inflow concentration of species B (Bin) to a non-zero value (Figure 8).

dyn("models/twoSpecies_basic.xlsx", update=c(h_B=0.36, Bin=1e5))

Figure 8: Prediction of the two-species model when the inferior competitor is subject to permanent invasion.

Invasion apparently allows for co-existence, however, we must be careful with wording. What we see in Figure 8 is not an example of stable co-existence. It is rather an example of co-presence!

In real world settings, it is actually quite difficult to distinguish cases of stable co-existence and co-presence because most ecosystems are interconnected and invasion is often a rule rather than an exception. Moreover, the possible time span of scientific investigations is usually too short to really test for stable co-existence. This is especially so when generation times are rather long and the differences in competitive ability are small.

4.3 Temporal variability can guard against extinction

Let’s study a scenario where each of the two competitors actually has its ecological niche. We assume that species A has a somewhat higher intrinsic growth rate at the reference temperature of 20°C compared to B. At the same time, we assume that the growth rate of B is more sensitive to temperature (Figure 9).

Tref <- 20                    # reference temperature
mu_Tref <- c(A=1, B=0.95)     # growth rates at reference temperature
theta <- c(A=1.05, B=1.1)     # temperature dependency parameters
                              # 1.07 is equivalent to factor 2 per 10 K

Tm <- 5:30                    # a gradient of temperatures

# temperature-adjusted growth rates of the competitors
mu_A <- mu_Tref["A"] * theta["A"]^(Tm - Tref)
mu_B <- mu_Tref["B"] * theta["B"]^(Tm - Tref)
plot(range(Tm), range(c(mu_A, mu_B)), type="n",
  xlab="Temperature", ylab="Temperature-adjusted growth rate")
lines(Tm, mu_A)
lines(Tm, mu_B, lty=2)
legend("topleft", bty="n", lty=1:2, legend=c("A","B"), title="Species")

Figure 9: Assumption on the temperature dependency of the growth rates of the competing species A and B. The particular mathematical formulation of temperature dependence is used by the famous WASP water quality model.

The extended model is specified by Table 5 and Table 6. Note the additional temperature-adjustment terms multiplied with the growth rates.

Table 5: State variables and parameters of the two species bioreactor model with contrasting temperature-dependent growth rates.

type	name	unit	description	default
variable	R	mg / mL	resource concentration	0
variable	A	cells / mL	abundance, organism A	100000
variable	B	cells / mL	abundance, organism B	100000
parameter	mu_A	1 / h	max. growth rate of A at ref. Temp.	1
parameter	mu_B	1 / h	max. growth rate of B at ref. Temp.	0.95
parameter	th_A	-	Temp. sensitivity parameter	1.05
parameter	th_B	-	Temp. sensitivity parameter	1.1
parameter	h_A	mg / mL	half saturation const. of A	0.18
parameter	h_B	mg / mL	half saturation const. of B	0.18
parameter	y_A	cells / mg	yield of A on R	1000000
parameter	y_B	cells / mg	yield of B on R	1000000
parameter	q	mL / h	flow rate	10
parameter	v	mL	reactor volume	100
parameter	Ain	cells / mL	abundance of A in input	0
parameter	Bin	cells / mL	abundance of B in input	0
parameter	Rin	mg / mL	input resource conc.	15
parameter	Tref	°C	Reference temp. for growth rate const.	20
function	Temp	°C	Time-varying system temperature	NA

Table 6: Simultaneous differential equations of the two species bioreactor model with contrasting temperature-dependent growth rates.

name	unit	description	rate	A	B	R
flow	1 / h	im-/export	q / v	Ain-A	Bin-B	Rin - R
growth_A	cells / mL / h	growth	mu_A * th_A^(Temp(time) - Tref) * R / (R + h_A) * A	1		-1 / y_A
growth_B	cells / mL / h	growth	mu_B * th_B^(Temp(time) - Tref) * R / (R + h_B) * B		1	-1 / y_B

Considering Figure 9, we expect A to win the competition game in a cold environment. This is confirmed by the following scenario (Figure 10).

Temp <- function(time) {17}                 # constant in time
dyn("models/twoSpecies_variability.xlsx",
  times=seq(0, 1200, 24))

Figure 10: Outcome of competition for cool conditions.

Likewise, we expect from Figure 9 that switching to a warm environment reverses the outcome of the competition experiment such that A goes extinct in the long run (Figure 11).

Temp <- function(time) {23}                 # constant in time
dyn("models/twoSpecies_variability.xlsx",
  times=seq(0, 1200, 24))

Figure 11: Outcome of competition for warm conditions.

Finally, we consider a scenario where the ambient temperature varies around the critical point where the two graphs of Figure 9 intersect. The simplest mean to achieve such fluctuations is by using periodic functions like \(sine\) or \(cosine\) as illustrated below (Figure 12).

Temp <- function(time, mini=17, maxi=25, lambda=72) {  # cyclic fluctuations
  pi <- 3.1415
  mini + (maxi - mini) * 0.5 * (1 + sin(time/lambda * 2*pi))
}
plot(0:120, Temp(0:120), type="l", xlab="Time", ylab="Temperature")
abline(v=seq(0, 240, 24), lty=3)

Figure 12: Temperature variation modeled by a sinus function. Range / amplitude was chosen to vary around the temperature of equal fitness from Figure 9. The frequency is determined by the `lambda` parameter in the corresponding code.

Now let’s rerun the competition model with the time-variable temperature forcing. To check for convergence against a steady state, we need to extend the time axis a bit (Figure 13). Note that calling the std() function is not a good idea now, because the state variables will never converge to constant values.

dyn("models/twoSpecies_variability.xlsx",
  times=seq(0, 24000, 480))

Figure 13: Outcome of competition for alternating temperatures (periodic forcing).

Under the particular settings chosen, stable coexistence seems possible. However, the range of suitable settings is quite narrow. As soon as one of the species profits from its preferred temperature over longer times, it is likely to outcompete the other one.

Nevertheless, temporal variability is certainly an important mechanism explaining species coexistence. As demonstrated by the example, it even works in very simply structured environments with no spatial structure and just a single resource where ecological niches can only be exploited in the time domain.

Even when temporal variability alone does not explain coexistence sufficiently, it can make an important contribution by delaying extinction at least.

4.4 Compensation of a fitness disadvange 1: Stabilizing effect of predation

Prediction the effects of predation is a surprisingly complex topic and you likely know the classical work of Lotka and Volterra. Here, we will focus on the simplest possible case where only one of the two species in the model reactor undergoes predation und thus experiences additional losses (Figure 14). We impose these additional loss on the species with otherwise greater fitness. The other (inferior) species is simply assumed to be inedible for any reason (may it be due to size, taste, swimming speed, …). The idea behind is that the two contrasting advantages / disadvantages may compensate each other.

Figure 14: Symbolic stoichiometry matrix of the two species bioreactor model with predation on a competitor.

The model is specified by Table 7 and Table 8. Note that the dependency of predation on prey density was considered by an extended Monod-term with an offset. The idea is that the predator becomes inactive when there is too few prey because energy gain does not compensate for energy spent on hunting or filtration activity.

Table 7: State variables and parameters of the two species bioreactor model with predation on a competitor.

type	name	unit	description	default
variable	R	mg / mL	resource concentration	0
variable	A	cells / mL	abundance, organism A	100000
variable	B	cells / mL	abundance, organism B	100000
variable	P	cells / mL	abundance of predator	100
parameter	mu_A	1 / h	max. growth rate of A	2
parameter	mu_B	1 / h	max. growth rate of B	1
parameter	mu_P	1 / h	max. growth rate of predator	0.1
parameter	h_A	mg / mL	half saturation const. of A	0.18
parameter	h_B	mg / mL	half saturation const. of B	0.18
parameter	h_P	cells / mL	half saturation const. of predator	1000000
parameter	y_A	cells / mg	yield of A on R	1000000
parameter	y_B	cells / mg	yield of B on R	1000000
parameter	y_P	cells / cells	yield of P on A or B	0.01
parameter	i_P	-	minimum prey level (relative to h_P)	0.01
parameter	q	mL / h	flow rate	2
parameter	v	mL	reactor volume	100
parameter	Ain	cells / mL	abundance of A in input	0
parameter	Bin	cells / mL	abundance of B in input	0
parameter	Pin	cells / mL	abundance of predator in input	0
parameter	Rin	mg / mL	input resource conc.	15
function	max	-	built-in function	NA

Table 8: Simultaneous differential equations of the two species bioreactor model with predation on a competitor.

name	unit	description	rate	A	B	P	R
flow	1 / h	im-/export	q / v	Ain-A	Bin-B	Pin-P	Rin - R
growth_A	cells / mL / h	growth	mu_A * R / (R + h_A) * A	1			-1 / y_A
growth_B	cells / mL / h	growth	mu_B * R / (R + h_B) * B		1		-1 / y_B
growth_P	cells / mL / h	growth	mu_P * max(0, (A - h_P * i_P) / (A - h_P * i_P + h_P)) * P	-1 / y_P		1

I already adjusted the default values of all parameters to demonstrate a stabilizing effect of predation on co-existence (Figure 15).

dyn("models/twoSpecies_predation.xlsx",
  times=seq(0, 4800, 40), maxcols=2)

Figure 15: Virtual competition experiment with predation pressure solely on the superior competitor.

For different parameter settings, you have to expect different outcomes, including extinction of one of the competitors. Finally, the example below demonstrates that predation does not necessarily trigger an infinite series of oscillations (Figure 16). It all depends on the relations between the various model parameters.

dyn("models/twoSpecies_predation.xlsx", update=c(mu_P=0.3),
  times=seq(0, 4800, 40), maxcols=2)

Figure 16: Virtual competition experiment with predation pressure solely on the superior competitor.

Finally, you should note that the oscillations visible in Figure 15 and Figure 16 are of a totally different kind compared to those seen in Figure 13. In the later figure, oscillations were due to the prescribed periodic temperature forcing imposed on the system. By contrast, oscillations in the predation example represent an internal feature of the system! The appear even if forcings are perfectly constants.

Now, you can probably imagine that the interpretation of oscillation in real-world observation data is challenging. In general, one has to expect that external variations at multiple time scales (daily, seasonal, other, yet other) possibly get superimposed by intrinsic oscillations generated by the system itself.

4.5 Compensation of a fitness disadvange 2: Persistence

In the previous example, the lower fitness of species B was compensated by increased losses of the competitor species A owing to its susceptibility to predation. The present example is somewhat similar but assumes an additional trait in species B rather than the presence of a predator. Specifically, it is assumed that species B tends to form biofilm covering the surface of the reactor (Figure 17, attached cells are denoted Ba). Thereby, species B partly escapes the permanent loss of suspended individuals from the reactor.

Figure 17: Symbolic stoichiometry matrix of the two species bioreactor model with biofilm formation.

The model is specified by Table 9 and Table 10. Note the two simultaneous process of attachment and detachment. In fact, the approach to model the distribution of species B over the suspended and biofilm phase is identical to what is used to describe chemical equilibria.

Table 9: State variables and parameters of the two species bioreactor model with biofilm formation.

type	name	unit	description	default
variable	R	mg / mL	resource concentration	0
variable	A	cells / mL	abundance, organism A	100000
variable	B	cells / mL	abundance, organism B	100000
variable	Ba	cells / cm2	abundance of attached B	0
parameter	mu_A	1 / h	max. growth rate of A	2
parameter	mu_B	1 / h	max. growth rate of B	2
parameter	h_A	mg / mL	half saturation const. of A	0.18
parameter	h_B	mg / mL	half saturation const. of B	0.18
parameter	y_A	cells / mg	yield of A on R	1000000
parameter	y_B	cells / mg	yield of B on R	1000000
parameter	ka	1 / h	attachment rate constant	1
parameter	kd	1 / h	detachment rate constant	1
parameter	ccap_Ba	cells / cm2	upper limit of biofilm density	100000000
parameter	q	mL / h	flow rate	10
parameter	v	mL	reactor volume	100
parameter	a	cm2	surface area for attachment	100
parameter	Rin	mg / mL	input resource conc.	15
parameter	Ain	cells / mL	abundance of A in input	0
parameter	Bin	cells / mL	abundance of B in input	0
function	max	-	built-in function	NA

Table 10: Simultaneous differential equations of the two species bioreactor model with biofilm formation.

name	unit	description	rate	A	B	Ba	R
flow	1 / h	im-/export	q / v	Ain-A	Bin-B		Rin - R
growth_A	cells / mL / h	growth	mu_A * R / (R + h_A) * A	1			-1 / y_A
growth_B	cells / mL / h	growth	mu_B * R / (R + h_B) * B		1		-1 / y_B
growth_Ba	cells / cm2 / h	growth	mu_B * R / (R + h_B) * max(0, 1 - Ba/ccap_Ba) * Ba			1	-1 / y_B * a / v
attach	cells / mL / h	attachment	ka * B		-1	v / a
detach	cells / cm2 / h	detachment	kd * Ba		a / v	-1

When the model is run with defaults settings, species B outcompetes A because the growth rates were set equal (Figure 18).

dyn("models/twoSpecies_biofilm.xlsx", maxcols=2)

Figure 18: Predicted dynamics with default values.

However, by lowering the growth rate of B, we can get to a state where the two species co-exist in a stable manner (Figure 19). Because of the additional biofilm trait, the system in general reacts a bit slower and the steady state is only reached after longer time.

dyn("models/twoSpecies_biofilm.xlsx", c(mu_B=1),
  times=seq(0, 4800, 24), maxcols=2)

Figure 19: Predicted dynamics with very low.

Again, the ability of the inferior competitor to attach to surfaces may facilitate coexistence, but it depends on the particular relation between parameter values. For this example, one could identify a feasible range of conditions, for example, by computing steady state solutions over a range of assumed growth rates of species B (Figure 20).

x <- NULL
for (mu_B in seq(0.7, 1.3, 0.1)) {
  x <- rbind(x, c(mu_B=mu_B,
    std("models/twoSpecies_biofilm.xlsx", c(mu_B=mu_B))
  ))
}
plot(range(x[,"mu_B"]), range(x[,c("A","B")]), type="n",
  xlab="mu_B", ylab="Steady state abundance (cells / ml)")
points(A ~ mu_B, data=x, pch="A")
points(B ~ mu_B, data=x, pch="B")

Figure 20: Predicted outcomes of competition for a range of growth rates of the inferior competitor. Letters denote the predicted steady state abundances for the respective types of suspended cells.

5 Summary

The mere presence of multiple species in a system at a certain point in time does not necessarily indicate their stable co-existence. Species presence may well be a result of invasion.
Modeling provides a unique means to study the conditions that allow for stable species co-existence.
In general, in accordance with theory, co-existence requires niches in the sense that each species is superior to competitors in a certain aspect.
Superiority arises from either faster growth or reduced losses (or both).
In real ecosystems, the co-presence of many seemingly similar species is often a result of multiple factors, including niches, invasion, spatial separation and mixing, temporary escape from competition, and temporal variability of boundary conditions.