Skip to content

Hyperbolic chemotaxis of cells and a chemical

Cells and a chemical moving with finite speed, whose two bumps of cells spread and attract each other.

  • Name: hyperbolic_chemotaxis
  • Field: biology/physiology
  • Formalism: ordinary differential equations
  • References: Bellomo et al. 2017

Select this system

from skcomplex.systems import create_simulation, create_system, describe

system = create_system('hyperbolic_chemotaxis')
entry = describe('hyperbolic_chemotaxis')
trajectory = create_simulation('hyperbolic_chemotaxis').prepare().run().trajectory_

Factory

@register_system(
    field='biology/physiology',
    tags=('chemotaxis', 'hyperbolic', 'kinetic theory', 'field', 'pde'),
    references=('Bellomo et al. 2017',),
    defaults={'until': 0.015, 'dt': 0.0005, 'solver': ScipyODESolver(method='BDF')},
)
def hyperbolic_chemotaxis(
    sites: int = 40,
    *,
    s: float = 100.0,
    D_n: float = 1.0,
    D_c: float = 0.001,
    alpha: float = 0.33,
    tau: float = 1.0,
    n0: float = 0.25,
    width: float = 0.03,
) -> System:
    """Return cells and a chemical moving with finite speed, whose two bumps of cells spread and attract each other.

    The model is the hyperbolic limit of a kinetic model of cells and a chemical that change velocity at random, with
    the speed `s` and turning rates set by the diffusions `D_n = s^2 / (2 mu1)` and `D_c = s^2 / (2 sigma)` of its
    parabolic limit, Keller-Segel, and `mu2 = mu1`. Following the test of Bellomo et al. (2017, section 6.6.2), the
    cells start as two Gaussian bumps of mass `n0` and width `width` at `(3 width, 3 width)` and `(-3 width, -3 width)`
    in the square `[-0.4, 0.4]^2`, here with periodic edges, the cells secrete the chemical, and all fluxes start at
    rest. The fluxes relax fast, so the bumps spread out at a finite speed while the cell mass `2 n0` is conserved and
    the chemical accumulates where the cells are.

    Args:
        sites:
            The number of sites along each side of the square.

        s:
            The speed of the cells and of the chemical.

        D_n:
            The diffusion of the cells in the parabolic limit.

        D_c:
            The diffusion of the chemical in the parabolic limit.

        alpha:
            The chemotactic sensitivity of the cells.

        tau:
            The time scale of the chemical.

        n0:
            The mass of each bump of cells.

        width:
            The width of each bump of cells.

    Returns:
        system:
            The system, with the cell density `tissue.n`, the chemical `tissue.c`, and their fluxes.

    Examples:
        >>> from skcomplex.systems.biology.physiology import hyperbolic_chemotaxis
        >>> round(hyperbolic_chemotaxis(20).elements[0].values['mu1'])
        5000
    """
    spacing = 0.8 / sites
    axis = -0.4 + spacing * (np.arange(sites) + 0.5)
    x, y = np.meshgrid(axis, axis, indexing='ij')
    offset = 3 * width
    bumps = np.exp(-((x - offset) ** 2 + (y - offset) ** 2) / (2 * width**2))
    bumps += np.exp(-((x + offset) ** 2 + (y + offset) ** 2) / (2 * width**2))
    turning = s**2 / (2 * D_n)
    values = {
        'n': n0 / (2 * np.pi * width**2) * bumps,
        's': s,
        'mu1': turning,
        'mu2': turning,
        'sigma': s**2 / (2 * D_c),
        'tau': tau,
        'alpha': alpha,
    }
    space = LatticeSpace((sites, sites), boundary='periodic', spacing=spacing)
    return System(elements=[FieldElement('tissue', element=ChemotacticTissueElement, space=space, values=values)])

Components

Elements:

  • ChemotacticTissueElement, from Element
class ChemotacticTissueElement(Element):
    """Site of a tissue where cells and a chemical move with finite speeds, the hyperbolic limit of velocity jumps.

    The cell density `n` and its flux `q`, and the chemical density `c` and its flux `j`, follow Eq. (6.50) of Bellomo
    et al. (2017) in two dimensions with the cells as the only source of the chemical, `dn/dt = -div q`, `dq/dt = -(s^2
    / 2) grad n - mu1 q + mu2 alpha n grad c`, `tau dc/dt = -div j + n`, and `tau dj/dt = -(s^2 / 2) grad c - sigma
    j`. Cells and chemical move at the speed `s` and turn at the rates `mu1` and `sigma`, and the cells turn towards
    the gradient of the chemical with the sensitivity `alpha`. In the limit of fast turning it is the Keller-Segel model
    with the diffusions `s^2 / (2 mu1)` and `s^2 / (2 sigma)`.

    Args:
        label:
            The label of the tissue.

        visual:
            How the element is drawn, or `None` for `default_visual_`.

        n:
            The initial cell density.

        qx:
            The initial flux of cells along the first axis.

        qy:
            The initial flux of cells along the second axis.

        c:
            The initial chemical density.

        jx:
            The initial flux of the chemical along the first axis.

        jy:
            The initial flux of the chemical along the second axis.

        s:
            The speed of the cells and of the chemical.

        mu1:
            The turning rate of the cells.

        mu2:
            The rate of turning of the cells towards the chemical.

        sigma:
            The turning rate of the chemical.

        tau:
            The time scale of the chemical.

        alpha:
            The chemotactic sensitivity of the cells.

    Examples:
        >>> from skcomplex.systems.biology.physiology import ChemotacticTissueElement
        >>> ChemotacticTissueElement('tissue').alpha
        0.33
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = None,
        *,
        n: Annotated[float, StateSpec()] = 0.0,
        qx: Annotated[float, StateSpec()] = 0.0,
        qy: Annotated[float, StateSpec()] = 0.0,
        c: Annotated[float, StateSpec()] = 0.0,
        jx: Annotated[float, StateSpec()] = 0.0,
        jy: Annotated[float, StateSpec()] = 0.0,
        s: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 100.0,
        mu1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 5000.0,
        mu2: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 5000.0,
        sigma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 5.0e6,
        tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
        alpha: Annotated[float, ParameterSpec()] = 0.33,
    ) -> None:
        super().__init__(label, visual=visual)
        self.n = n
        self.qx = qx
        self.qy = qy
        self.c = c
        self.jx = jx
        self.jy = jy
        self.s = s
        self.mu1 = mu1
        self.mu2 = mu2
        self.sigma = sigma
        self.tau = tau
        self.alpha = alpha

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the balance laws of the cells and the chemical with their turning relaxations.

        Args:
            s:
                The namespace of the symbols the laws can use.

        Returns:
            laws:
                The laws of the component.
        """
        pressure = s.s**2 / 2
        cells, chemical = gradient(s.n), gradient(s.c)
        return [
            RateLaw(s.n, -divergence([s.qx, s.qy])),
            RateLaw(s.qx, -pressure * cells[0] - s.mu1 * s.qx + s.mu2 * s.alpha * s.n * chemical[0]),
            RateLaw(s.qy, -pressure * cells[1] - s.mu1 * s.qy + s.mu2 * s.alpha * s.n * chemical[1]),
            RateLaw(s.c, (-divergence([s.jx, s.jy]) + s.n) / s.tau),
            RateLaw(s.jx, (-pressure * chemical[0] - s.sigma * s.jx) / s.tau),
            RateLaw(s.jy, (-pressure * chemical[1] - s.sigma * s.jy) / s.tau),
        ]

Equations

tissue.n[i]: diff(n_tissue(t), t) == -partial(qx_tissue(t), 0) - partial(qy_tissue(t), 1)
tissue.qx[i]: diff(qx_tissue(t), t) == alpha_tissue*mu2_tissue*n_tissue(t)*partial(c_tissue(t), 0) - 1/2*s_tissue^2*partial(n_tissue(t), 0) - mu1_tissue*qx_tissue(t)
tissue.qy[i]: diff(qy_tissue(t), t) == alpha_tissue*mu2_tissue*n_tissue(t)*partial(c_tissue(t), 1) - 1/2*s_tissue^2*partial(n_tissue(t), 1) - mu1_tissue*qy_tissue(t)
tissue.c[i]: diff(c_tissue(t), t) == (n_tissue(t) - partial(jx_tissue(t), 0) - partial(jy_tissue(t), 1))/tau_tissue
tissue.jx[i]: diff(jx_tissue(t), t) == -1/2*(s_tissue^2*partial(c_tissue(t), 0) + 2*sigma_tissue*jx_tissue(t))/tau_tissue
tissue.jy[i]: diff(jy_tissue(t), t) == -1/2*(s_tissue^2*partial(c_tissue(t), 1) + 2*sigma_tissue*jy_tissue(t))/tau_tissue