Immune competition
The kinetic model of the competition between tumour cells and immune cells.
- Name:
immune_competition - Field:
biology/oncology - Formalism: ordinary differential equations
- References: Bellomo et al. 2017; Bellomo and Delitala 2008; Bellomo, Bellouquid, and Delitala 2004
Factory
@register_model(
field='biology/oncology',
tags=('tumour', 'immune', 'kinetic', 'distribution'),
references=IMMUNE_REFERENCES,
defaults={'until': 40, 'dt': 0.5},
observables=lambda system: immune_observables(_classes_of(system)),
)
def immune_competition(
classes: int = 5,
*,
activation: float = 0.75,
tumour: float = 0.05,
immune: float = 1.0,
growth: float = 0.3,
progression: float = 0.05,
kill: float = 1.0,
inhibition: float = 0.2,
) -> System:
"""Return the kinetic model of the competition between tumour cells and immune cells.
Tumour cells and immune cells are distributions over activity classes, the progression of tumour cells and the
activation of immune cells, held by site groups on a lattice of the classes. Tumour cells proliferate in the host
and progress, active immune cells destroy them, and tumour cells inhibit the activation of immune cells. The
tumour starts in the lowest progression class and the immune cells in the class nearest `activation`. While
`kill activation immune` exceeds `growth`, the immune response wins and the tumour density decays. Below that
threshold the tumour grows, inhibits the immune cells further, and escapes.
Args:
classes:
The number of activity classes of both distributions.
activation:
The initial activation of the immune cells, between 0 and 1.
tumour:
The initial density of tumour cells.
immune:
The density of immune cells.
growth:
The proliferation rate of tumour cells.
progression:
The rate at which tumour cells progress one class.
kill:
The rate at which fully active immune cells destroy tumour cells.
inhibition:
The rate at which the most progressed tumour cells inhibit immune cells.
Returns:
system:
The model, with the distributions `tumour` and `immune`, the host `host`, and the interactions `growth`
and `competition`. Its standard observables are the tumour density, the mean tumour progression, and the
mean immune activation.
Examples:
>>> from skcomplex.models.biology.oncology import immune_competition
>>> [item.label for item in immune_competition().elements]
['tumour', 'immune', 'host']
"""
activities = activity_grid(classes)
activation_rates = np.tile(activities, (classes, 1))
tumour_density = np.zeros(classes)
tumour_density[0] = tumour
immune_density = np.zeros(classes)
immune_density[int(np.argmin(np.abs(activities - activation)))] = immune
return System(
elements=[
GroupElement('tumour', element=DensityElement, size=classes, on='sites', values={'f': tumour_density}),
GroupElement('immune', element=DensityElement, size=classes, on='sites', values={'f': immune_density}),
GroupElement('host', element=DensityElement, size=1, values={'f': 1.0}),
],
interactions=[
EncounterInteraction(
'proliferation',
candidate='tumour',
field='host',
rate=lambda _: growth,
outcome='proliferative',
),
EncounterInteraction(
'progression',
candidate='tumour',
field='host',
rate=lambda _: progression,
kernel=lambda _: step_table(classes, 1),
),
EncounterInteraction(
'destruction',
candidate='tumour',
field='immune',
rate=lambda _: activation_rates * kill,
outcome='destructive',
),
EncounterInteraction(
'inhibition',
candidate='immune',
field='tumour',
rate=lambda _: activation_rates * inhibition,
kernel=lambda _: step_table(classes, -1, classes),
),
],
space=LatticeSpace(classes, boundary='fixed'),
time=ContinuousTime(),
)
Components
class DensityElement(Element):
"""Value of a density at one class, a member of a distribution over the classes.
A distribution is a group of densities on the sites of a lattice of classes, one member per class, and encounter,
proliferation, and destruction laws change it. A density has no laws of its own.
Args:
label:
The label of the density, unique among the parts of its system.
visual:
How the element is drawn, or `None` for `default_visual_`.
f:
The initial density, a nonnegative number.
Attributes:
default_visual_:
The visual of the element when it sets none.
Examples:
>>> from skcomplex.core.elements import DensityElement, GroupElement
>>> from skcomplex.core.spaces import LatticeSpace
>>> classes = LatticeSpace(5, boundary='fixed')
>>> group = GroupElement('distribution', element=DensityElement, space=classes, on='sites', values={'f': 0.2})
>>> group.check().prepare().size_
5
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
f: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
super().__init__(label, visual=visual)
self.f = f
class EncounterInteraction(Interaction):
"""Encounters of a candidate density with a field density over activity classes, which move, remove, or create.
The candidate and the field are groups of `DensityElement` on the sites of lattices of activity classes, one member
per class. A candidate at class `h` meets a field at class `k` at the rate `eta[h, k]`. In the outcome
`'conservative'` the candidate moves to class `i` with the probability `kernel[h, k, i]`, in the outcome
`'destructive'` it is removed, and in the outcome `'proliferative'` it creates a new candidate at class `i` with
the probability `kernel[h, k, i]`, or at its own class without a kernel. The rate and the kernel are functions of
the namespace, which holds `namespace.candidate`, `namespace.field`, and the parameters of the interaction.
On a `PhaseSpace` the classes are the points of the class axes, and encounters pair classes at the same position. A
rate that is one expression may read the `coordinate` of the class of the candidate. Without a field the encounter
is a linear redistribution: a conservative candidate at class `h` moves to class `i` at the rate `eta[h]` with the
probability `kernel[h, i]`.
With a `cross_section` the candidate is a field on a `PhaseSpace` with two interval class axes, and its rate at
every position is the binary collision integral over those axes, evaluated by the fast spectral method on a padded
grid.
The density is conserved exactly, and the momentum and the energy within the spectral accuracy of the method.
Args:
label:
The label of the interaction, unique among the parts of its system.
candidate:
The path of the group of the candidate density.
field:
The path of the group of the field density, or `None` for a linear redistribution.
rate:
The function `rate(namespace)` of the encounter rate, a number or an expression, or an array `eta[h, k]`
of them, or `eta[h]` without a field, or `None` with a cross-section.
kernel:
The function `kernel(namespace)` of the array `kernel[h, k, i]` of the probabilities of the class of the
moved or created candidate, or `kernel[h, i]` without a field, or `None` for a destructive encounter or the
candidate's own class.
outcome:
`'conservative'`, `'destructive'`, or `'proliferative'`.
state:
The name of the state of the densities.
cross_section:
The function `B(speed, cosine)` of arrays of the relative speed and the cosine of the deflection angle, or
`None` for encounters at a rate.
factor:
A function of the namespace that returns a factor of the collision frequency, such as a function of the
density moment, or `None` for none.
Raises:
ValidationError:
The outcome is unknown, the state is no name, a function is not callable, a conservative encounter has no
kernel, an encounter without a field is not conservative, an encounter has neither a rate nor a
cross-section, or a cross-section comes with a rate, a kernel, a field, or an outcome.
Examples:
>>> import numpy as np
>>> from skcomplex.core.interactions import EncounterInteraction
>>> kernel = np.full((3, 3, 3), 1 / 3)
>>> mixing = EncounterInteraction(
... 'mixing', candidate='population', field='population', rate=lambda s: 1.0, kernel=lambda s: kernel
... )
>>> mixing.check().outcome
'conservative'
>>> collide = EncounterInteraction('collide', candidate='gas', cross_section=lambda speed, cosine: 1.0)
>>> collide.check().rate is None
True
"""
def __init__(
self,
label: str,
*,
candidate: Annotated[str, RoleSpec()],
field: Annotated[str | None, RoleSpec(optional=True)] = None,
rate: Callable[[SimpleNamespace], object] | None = None,
kernel: Callable[[SimpleNamespace], object] | None = None,
outcome: str = 'conservative',
state: str = 'f',
cross_section: Callable[..., object] | None = None,
factor: Callable[[SimpleNamespace], object] | None = None,
) -> None:
super().__init__(label)
self.candidate = candidate
self.field = field
self.rate = rate
self.kernel = kernel
self.outcome = outcome
self.state = state
self.cross_section = cross_section
self.factor = factor
def _check_collision(self) -> None:
"""Raise an error for a cross-section or a factor that is no function, or a cross-section with others."""
if not callable(self.cross_section):
error_msg = f'The cross-section of {self.label!r} is {self.cross_section!r}. Give a function of two arrays.'
raise ValidationError(error_msg)
if self.factor is not None and not callable(self.factor):
error_msg = f'The factor of {self.label!r} is {self.factor!r}. Give a function of the namespace.'
raise ValidationError(error_msg)
others = self.rate, self.kernel, self.field
if any(item is not None for item in others) or self.outcome != 'conservative':
error_msg = f'The cross-section of {self.label!r} sets the rate. Give no rate, kernel, field, or outcome.'
raise ValidationError(error_msg)
def _check(self) -> None:
"""Raise an error for an unknown outcome, a function that is not callable, or a conservative kernel missing.
Raises:
ValidationError:
The outcome is unknown, the state is no name, a function is not callable, a conservative encounter has
no kernel, an encounter without a field is not conservative, an encounter has neither a rate nor a
cross-section, or a cross-section comes with a rate, a kernel, a field, or an outcome.
"""
super()._check()
if not isinstance(self.state, str) or not self.state.isidentifier():
error_msg = f'{type(self).__name__} {self.label!r} reads the state {self.state!r}, which is no name.'
raise ValidationError(error_msg)
if self.cross_section is not None:
self._check_collision()
return
if self.factor is not None:
error_msg = f'The factor of {self.label!r} scales a cross-section. Give `cross_section`.'
raise ValidationError(error_msg)
if self.outcome not in OUTCOMES:
error_msg = (
f'The outcome of encounter {self.label!r} is {self.outcome!r}. Use one of {", ".join(OUTCOMES)}.'
)
raise ValidationError(error_msg)
if not callable(self.rate) or (self.kernel is not None and not callable(self.kernel)):
error_msg = f'The rate and the kernel of encounter {self.label!r} must be functions of the namespace.'
raise ValidationError(error_msg)
if self.outcome == 'conservative' and self.kernel is None:
error_msg = f'The conservative encounter {self.label!r} moves candidates by a kernel. Give `kernel`.'
raise ValidationError(error_msg)
if self.field is None and self.outcome != 'conservative':
error_msg = f'The {self.outcome} encounter {self.label!r} reads a density of others. Give `field`.'
raise ValidationError(error_msg)
def _linear_laws(self, candidate: Expression, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the encounter law of a redistribution without a field, one field class of value one."""
assert self.rate is not None
rate = np.asarray(self.rate(namespace), dtype=object)
rates = rate if rate.ndim == 0 else rate.reshape(-1, 1)
kernel = np.asarray(self.kernel(namespace) if self.kernel is not None else (), dtype=object)
return [EncounterLaw(candidate, None, rates[()], np.expand_dims(kernel, axis=1), label=self.label)]
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the encounter, bilinear, or collision law of the candidate density.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.candidate` and `namespace.field`.
Returns:
laws:
The law of the encounters.
"""
candidate = getattr(self._get_member(namespace, 'candidate'), self.state)
if self.cross_section is not None:
factor = SR(1 if self.factor is None else self.factor(namespace))
return [
RateLaw(candidate, factor * collision(candidate), cross_section=self.cross_section, label=self.label),
]
if self.field is None:
return self._linear_laws(candidate, namespace)
field = getattr(self._get_member(namespace, 'field'), self.state)
assert self.rate is not None
rate = self.rate(namespace)
kernel = None if self.kernel is None else self.kernel(namespace)
if self.outcome == 'destructive':
loss = -np.asarray(rate, dtype=object)
return [BilinearLaw(candidate, candidate, field, loss, label=self.label)]
if self.outcome == 'proliferative':
return [BilinearLaw(candidate, candidate, field, rate, kernel, label=self.label)]
return [EncounterLaw(candidate, field, rate, kernel, label=self.label)]
Equations
tumour.f[i]: diff(f_tumour(t), t) == 0
immune.f[i]: diff(f_immune(t), t) == 0
host.f[i]: diff(f_host(t), t) == 0