Consumer-resource models of substitutable and essential resources
Consumer species that compete only through the resources they take up, by MacArthur's niche theory.
- Name:
consumer_resource - Field:
biology/ecology - Formalism: ordinary differential equations
- References: MacArthur 1970; Tilman 1982; Mehta and Blumenthal 2024; Gros 2015
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('consumer_resource')
entry = describe('consumer_resource')
trajectory = create_simulation('consumer_resource').prepare().run().trajectory_
Factory
@register_system(
field='biology/ecology',
tags=('consumer-resource', 'competitive exclusion', 'niche theory', 'ode'),
references=('MacArthur 1970', 'Tilman 1982', 'Mehta and Blumenthal 2024', 'Gros 2015'),
defaults={'until': 100.0, 'dt': 0.5},
)
def consumer_resource(
consumers: Sequence[float] | None = None,
resources: Sequence[float] | None = None,
*,
c: Sequence[Sequence[float]] = UPTAKE,
m: Sequence[float] | float = 1.0,
e: Sequence[float] | float = 1.0,
w: Sequence[float] | float = 1.0,
r: Sequence[float] | float = 1.0,
K: Sequence[float] | float = CAPACITIES,
kappa: Sequence[float] | float = CAPACITIES,
omega: Sequence[float] | float = 1.0,
supply: str = 'logistic',
kind: str = 'substitutable',
n: float = -10.0,
v: Sequence[Sequence[float]] | None = None,
) -> System:
"""Return consumer species that compete only through the resources they take up, by MacArthur's niche theory.
The `'logistic'` supply of substitutable resources is MacArthur's consumer-resource model, `dN_i/dt = e_i N_i
(sum_b w_b c_ib R_b - m_i)` and `dR_a/dt = r_a R_a (K_a - R_a) / K_a - sum_j c_ja N_j R_a`. The `'chemostat'` supply
`kappa_a - omega_a R_a` models nutrients that flow in from outside, and essential resources limit growth jointly
through `(sum_b (c_ib R_b)^n)^(1/n)`, consumed by the stoichiometry `v`. At most as many consumers as resources
coexist, at the point where their zero net growth isoclines `g_i(R) = m_i` cross. With the defaults, the settings of
Mehta and Blumenthal (2024, figure 4), both consumers coexist at `R = (5/3, 5/9)`, the point of the region where
no consumer can invade that is closest to the supply point `K`.
The competition for one resource of Gros (2015, exercise 8), `dx_i/dt = (f r_i - d) x_i` and `df/dt = a - f sum_i
r_i x_i`, is the chemostat with one resource, `c = r`, `m = d`, `kappa = a`, and `omega = 0`. The total population
tends to `a / d`, the resource to `min_i d / r_i`, and only the consumer with the largest `r_i` survives.
Args:
consumers:
The initial population sizes, or `None` for 0.1 each.
resources:
The initial resource abundances, or `None` for the capacities `K` or the inflows `kappa`.
c:
The uptake rates, one row per consumer and one column per resource.
m:
The maintenance energies of the consumers.
e:
The efficiencies of the consumers.
w:
The energies of the resources.
r:
The growth rates of self-renewing resources.
K:
The carrying capacities of self-renewing resources.
kappa:
The inflows of supplied resources.
omega:
The dilution rates of supplied resources.
supply:
`'logistic'` for self-renewing resources or `'chemostat'` for supplied ones.
kind:
`'substitutable'` or `'essential'` resources.
n:
The exponent of the growth on essential resources, below -1.
v:
The stoichiometry of essential resources, one row per consumer, or `None` for `c`.
Returns:
system:
The system, with the populations `consumers.N` and the resources `resources.R`.
Raises:
ValidationError:
The uptake rates are not one row per consumer with one column per resource.
Examples:
>>> from skcomplex.systems.biology.ecology import consumer_resource
>>> [group.prepare().size_ for group in consumer_resource().elements]
[2, 2]
"""
uptake = np.atleast_2d(np.asarray(c, dtype=float))
n_consumers, n_resources = uptake.shape
sizes = np.full(n_consumers, 0.1) if consumers is None else np.asarray(consumers, dtype=float)
supplied = K if supply == 'logistic' else kappa
abundances = np.asarray(resources if resources is not None else np.broadcast_to(supplied, n_resources), float)
if sizes.shape != (n_consumers,) or abundances.shape != (n_resources,):
error_msg = (
f'The uptake rates of shape {uptake.shape} need {n_consumers} initial populations and {n_resources} '
'initial resources.'
)
raise ValidationError(error_msg)
consumer_group = GroupElement(
'consumers',
element=ConsumerElement,
size=n_consumers,
values={'N': sizes, 'm': m, 'e': e},
)
resource_group = GroupElement(
'resources',
element=ResourceElement,
size=n_resources,
values={'R': np.array(abundances, dtype=float), 'w': w, 'r': r, 'K': K, 'kappa': kappa, 'omega': omega},
)
uptake_rates = tuple(map(tuple, uptake.tolist()))
growth = ConsumerGrowthInteraction(
'growth',
consumers='consumers',
resources='resources',
c=uptake_rates,
kind=kind,
n=n,
)
turnover = ResourceTurnoverInteraction(
'turnover',
consumers='consumers',
resources='resources',
c=uptake_rates,
kind=kind,
n=n,
supply=supply,
v=v,
)
return System(elements=[consumer_group, resource_group], interactions=[growth, turnover])
Components
Elements:
ConsumerElement, fromElementResourceElement, fromElement
Interactions:
ConsumerGrowthInteraction, fromInteractionResourceTurnoverInteraction, fromInteraction
class ConsumerElement(Element):
"""Consumer species that grows on the energy it takes from resources and needs a maintenance energy to persist.
Args:
label:
The label of the species or of the group of species.
visual:
How the element is drawn, or `None` for `default_visual_`.
N:
The initial population size.
e:
The efficiency that converts the energy taken in into growth.
m:
The maintenance energy, the intake at which the population neither grows nor declines.
Examples:
>>> from skcomplex.systems.biology.ecology import ConsumerElement
>>> ConsumerElement('consumers', m=0.5).m
0.5
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
N: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
e: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
m: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.N = N
self.e = e
self.m = m
class ResourceElement(Element):
"""Resource that consumers take up, either self-renewing with logistic growth or supplied from outside.
A self-renewing resource grows as `r R (K - R) / K`, and a supplied one flows in at `kappa` and is washed out at
the rate `omega`. The interaction that couples the resource to its consumers chooses which supply applies.
Args:
label:
The label of the resource or of the group of resources.
visual:
How the element is drawn, or `None` for `default_visual_`.
R:
The initial abundance.
w:
The energy a unit of the resource carries.
r:
The growth rate of a self-renewing resource at low abundance.
K:
The carrying capacity of a self-renewing resource.
kappa:
The inflow of a supplied resource.
omega:
The dilution rate of a supplied resource.
Examples:
>>> from skcomplex.systems.biology.ecology import ResourceElement
>>> ResourceElement('resources', K=4.8).K
4.8
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
R: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
w: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
r: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
K: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
kappa: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
omega: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.R = R
self.w = w
self.r = r
self.K = K
self.kappa = kappa
self.omega = omega
class ConsumerGrowthInteraction(ConsumptionInteraction):
"""Growth of consumers on the resources they take up, minus their maintenance.
On substitutable resources, `dN_i/dt = e_i N_i (sum_b w_b c_ib R_b - m_i)`. On essential resources, `dN_i/dt =
N_i (g_i - m_i)` with `g_i = (sum_b (c_ib R_b)^n)^(1/n)`.
Args:
label:
The label of the interaction.
consumers:
The path of the group of consumers.
resources:
The path of the group of resources.
c:
The uptake rates, one row per consumer and one column per resource.
kind:
`'substitutable'` or `'essential'` resources.
n:
The exponent of the growth on essential resources, below -1.
Examples:
>>> from skcomplex.systems.biology.ecology import ConsumerGrowthInteraction
>>> ConsumerGrowthInteraction('growth', consumers='consumers', resources='resources').c
((1.0,),)
"""
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the growth of the consumers.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
consumers, resources = s.consumers, s.resources
if self.kind == 'essential':
growth = self._essential_growth(resources.R) - consumers.m
else:
intake = total(resources.w * resources.R, weights=np.asarray(self.c, dtype=float))
growth = consumers.e * (intake - consumers.m)
return [RateLaw(consumers.N, consumers.N * growth)]
class ResourceTurnoverInteraction(ConsumptionInteraction):
"""Supply of resources, by logistic renewal or a chemostat, and their depletion by the consumers.
A `'logistic'` supply renews resource `a` at `r_a R_a (K_a - R_a) / K_a`, and a `'chemostat'` supply at `kappa_a -
omega_a R_a`. Substitutable resources lose `sum_j c_ja N_j R_a`, and essential resources lose `sum_j v_ja N_j g_j`
by the fixed stoichiometry `v` of the consumers.
Args:
label:
The label of the interaction.
consumers:
The path of the group of consumers.
resources:
The path of the group of resources.
c:
The uptake rates, one row per consumer and one column per resource.
kind:
`'substitutable'` or `'essential'` resources.
n:
The exponent of the growth on essential resources, below -1.
supply:
`'logistic'` for self-renewing resources or `'chemostat'` for resources supplied from outside.
v:
The stoichiometry of essential resources in the consumers, one row per consumer, or `None` for `c`.
Raises:
ValidationError:
The supply or the kind is unknown, or the exponent of essential resources is not below -1.
Examples:
>>> from skcomplex.systems.biology.ecology import ResourceTurnoverInteraction
>>> ResourceTurnoverInteraction('turnover', consumers='consumers', resources='resources').supply
'logistic'
"""
def __init__(
self,
label: str,
*,
consumers: Annotated[str, RoleSpec()],
resources: Annotated[str, RoleSpec()],
c: Sequence[Sequence[float]] = ((1.0,),),
kind: str = 'substitutable',
n: float = -10.0,
supply: str = 'logistic',
v: Sequence[Sequence[float]] | None = None,
) -> None:
super().__init__(label, consumers=consumers, resources=resources, c=c, kind=kind, n=n)
self.supply = supply
self.v = v
def _check(self) -> None:
"""Check the supply, the kind, and the exponent.
Raises:
ValidationError:
The supply or the kind is unknown, or the exponent of essential resources is not below -1.
"""
super()._check()
if self.supply not in SUPPLIES:
error_msg = f'The supply must be one of {", ".join(SUPPLIES)}, not {self.supply!r}.'
raise ValidationError(error_msg)
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the supply and the depletion of the resources.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
consumers, resources = s.consumers, s.resources
R = resources.R
if self.supply == 'logistic':
renewal = resources.r * R * (resources.K - R) / resources.K
else:
renewal = resources.kappa - resources.omega * R
if self.kind == 'essential':
stoichiometry = np.asarray(self.c if self.v is None else self.v, dtype=float)
depletion = total(consumers.N * self._essential_growth(R), weights=stoichiometry.T)
else:
depletion = R * total(consumers.N, weights=np.asarray(self.c, dtype=float).T)
return [RateLaw(R, renewal - depletion)]
Equations
consumers.N[i]: diff(N_consumers(t), t) == (_reduced_0 - m_consumers)*e_consumers*N_consumers(t)
resources.R[i]: diff(R_resources(t), t) == -_reduced_1*R_resources(t) + (K_resources - R_resources(t))*r_resources*R_resources(t)/K_resources