One-dimensional active viscoelastic and poroelastic gel with contractile cells
A one-dimensional active viscoelastic gel, where contractile cells concentrate by dragging the gel.
- Name:
active_viscoelastic_gel - Field:
physics/active_matter - Formalism: ordinary differential equations
- References: Brauns and Marchetti 2024; arXiv:2306.08868
Factory
@register_model(
field='physics/active_matter',
tags=('active gel', 'contractile', 'viscoelastic', 'field', 'inversion', 'pattern formation', 'pde'),
references=('Brauns and Marchetti 2024', 'arXiv:2306.08868'),
defaults={'until': 30.0, 'dt': 0.1, 'solver': RK4Solver(dt=0.005)},
observables=lambda system: [ObservableSpec('v', _gel_velocity, group='gel', reduce=_sites)],
)
def active_viscoelastic_gel(
n_sites: int = 128,
*,
length: float = 10.0,
density: float = 1.0,
perturbation: float = 0.01,
D: float = 0.1,
gamma: float = 1.0,
eta: float = 1.0,
E: float = 0.1,
T0: float = 4.0,
c_sat: float = 1.0,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return a one-dimensional active viscoelastic gel, where contractile cells concentrate by dragging the gel.
On a periodic line, the density `c` of the cells and the displacement `u` of a Kelvin-Voigt gel follow
c_t + (v c)_x = D c_xx, u_t = v, (gamma - eta d_xx) v = d_x (E u_x + Ta(c)),
with the active stress `Ta(c) = T0 c / (c + c_sat)`. The force balance holds at every step, as the modified
Helmholtz inversion of the stress gradient with the screening `sqrt(gamma / eta)`. Without stiffness the uniform
state loses stability at long wavelengths above the contractile threshold `c Ta'(c) > gamma D`, and the density
concentrates into peaks. The source leaves the active stress open, and the saturating one keeps the peaks finite.
The catalog entry streams the velocity as the observable `v`.
Args:
n_sites:
The number of sites.
length:
The length of the line.
density:
The mean density of the cells.
perturbation:
The standard deviation of the initial density around its mean.
D:
The diffusion coefficient of the cells.
gamma:
The friction of the gel with the substrate.
eta:
The viscosity of the gel.
E:
The stiffness of the gel.
T0:
The largest active stress of the cells.
c_sat:
The density at which the active stress is half its largest value.
random_state:
The seed or generator of the initial density.
Returns:
system:
The model, with the density `gel.c` and the displacement `gel.u`.
Examples:
>>> from skcomplex.models.physics.active_matter import active_viscoelastic_gel
>>> active_viscoelastic_gel(16).elements[0].values['eta']
1.0
"""
rng = np.random.default_rng(random_state)
values = {
'c': density + perturbation * rng.standard_normal(n_sites),
'u': 0.0,
'D': D,
'gamma': gamma,
'eta': eta,
'E': E,
'T0': T0,
'c_sat': c_sat,
}
gel = FieldElement(
'gel',
element=ContractileGelSiteElement,
space=LatticeSpace(n_sites, spacing=length / n_sites),
values=values,
)
transport = AdvectionInteraction('transport', target='gel', state='c', drift=_transport, limiter=None, riemann=None)
diffusion = DiffusionInteraction('diffusion', target='gel', state='c', coefficient=_diffusion)
displacement = RateInteraction('displacement', target='gel', state='u', rate=_transport)
return System(elements=[gel], interactions=[transport, diffusion, displacement])
Components
class ContractileGelSiteElement(Element):
"""Site of a one-dimensional active viscoelastic gel, the density of contractile cells and the gel displacement.
The site holds the density `c` of the cells and the displacement `u` of a Kelvin-Voigt gel of stiffness `E` and
viscosity `eta`, with friction `gamma` against a rigid substrate. The cells diffuse with the coefficient `D` and
exert the active stress `T0 c / (c + c_sat)`.
Args:
label:
The label of the site or the gel.
visual:
How the site is drawn, or `None` for `default_visual_`.
c:
The initial density of the cells.
u:
The initial displacement of the gel.
D:
The diffusion coefficient of the cells.
gamma:
The friction of the gel with the substrate.
eta:
The viscosity of the gel.
E:
The stiffness of the gel.
T0:
The largest active stress of the cells.
c_sat:
The density at which the active stress is half its largest value.
Examples:
>>> from skcomplex.models.physics.active_matter import ContractileGelSiteElement
>>> ContractileGelSiteElement('site').T0
4.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
c: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
u: Annotated[float, StateSpec()] = 0.0,
D: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.1,
gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
eta: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
E: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
T0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 4.0,
c_sat: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.c = c
self.u = u
self.D = D
self.gamma = gamma
self.eta = eta
self.E = E
self.T0 = T0
self.c_sat = c_sat
class AdvectionInteraction(ConservationInteraction):
"""Advection `du/dt + divergence(a u) = source(ns)` of a state `u` of a field by a drift field `a(ns)`.
The flux of the conservation law is the drift times the state, written in conservative form, so the total of the
state on a periodic field stays constant and, on a spectral field, the dealiasing of products applies to it. The
drift is one expression for a field of one dimension, or one per dimension, of the states and parameters of the
field, and may read other states than the one it carries, such as the flow of a fluid derived from a
streamfunction, or the drift of a probability density in a Fokker-Planck equation.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the field whose state is carried.
state:
The name of the carried state of the field.
drift:
A function of the namespace that returns the drift, or one drift per dimension.
source:
A function of the namespace that returns the source, or `None` for none.
limiter:
The slope limiter of the reconstruction of a finite-volume field, `'minmod'`, `'van-leer'`, or `None`.
riemann:
The numerical flux of a finite-volume field, `'rusanov'` or `'godunov'`, or `None` for the divergence of
the flux on a field of another discretization.
sources:
The paths of other elements the drift and the source read, as the namespaces `sources`. A field over the
position lattice of a target on a phase space is read broadcast along the class axes.
Raises:
ValidationError:
The state is not a name, a function is not callable, or the limiter or the numerical flux is unknown.
Examples:
>>> from typing import Annotated
>>> import numpy as np
>>> from skcomplex.core.elements import Element, FieldElement
>>> from skcomplex.core.interactions import AdvectionInteraction
>>> from skcomplex.core.spaces import LatticeSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class CellElement(Element):
... def __init__(self, label, u: Annotated[float, StateSpec()] = 0.0):
... super().__init__(label)
... self.u = u
>>> space = LatticeSpace((50,), spacing=0.1, boundary='periodic')
>>> values = {'u': np.exp(-((np.arange(50) - 25) / 5.0) ** 2)}
>>> volumes = 'finite-volume'
>>> field = FieldElement('tracer', element=CellElement, space=space, values=values, discretization=volumes)
>>> drift = AdvectionInteraction('drift', target='tracer', state='u', drift=lambda s: 1.0)
>>> system = System(elements=[field], interactions=[drift])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=1.0).trajectory_
>>> total = trajectory.values('tracer.u')[-1].sum()
>>> round(float(total), 6) == round(float(values['u'].sum()), 6)
True
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'u',
drift: Callable[[SimpleNamespace], object] | None = None,
source: Callable[[SimpleNamespace], object] | None = None,
limiter: str | None = 'minmod',
riemann: str | None = 'rusanov',
sources: Annotated[Sequence[str], RoleSpec(variadic=True)] = (),
) -> None:
super().__init__(
label,
target=target,
state=state,
source=source,
limiter=limiter,
riemann=riemann,
sources=sources,
)
self.drift = drift
def _check(self) -> None:
"""Raise an error for a state that is not a name, a drift that is not callable, or an unknown scheme.
Raises:
ValidationError:
The type declares states, the state is not a name, a function is not callable, or the limiter or the
numerical flux is unknown.
"""
super()._check()
if not callable(self.drift):
error_msg = f'The drift of {self.label!r} is {self.drift!r}. Give a function of the namespace.'
raise ValidationError(error_msg)
def _get_flux(self, namespace: SimpleNamespace) -> object:
"""Return the flux `a u`, one component per component of the drift.
Args:
namespace:
The namespace of the symbols the flux can use, with the namespace of the field.
Returns:
flux:
The flux, one expression or one per dimension.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
drift = self.drift(namespace) if callable(self.drift) else None
if isinstance(drift, FreeModuleElement | list | tuple):
return [SR(component) * state for component in drift]
return SR(drift) * state
class DiffusionInteraction(ConservationInteraction):
"""Diffusion `du/dt = D laplacian(u)` of a state `u` of a field, the conservation law of the flux `-D gradient(u)`.
The coefficient `D` is uniform in space, a number or a function of the namespace that returns an expression of
parameters and time. The interaction lowers to the Laplacian of the field, which its discretization evaluates, so
the total of the state on a periodic field stays constant and a point mass spreads as the heat kernel.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the field whose state diffuses.
state:
The name of the state of the field that diffuses.
coefficient:
A function of the namespace that returns the diffusion coefficient, or a number.
Raises:
ValidationError:
The state is not a name, or the coefficient is neither callable nor a number.
Examples:
>>> from typing import Annotated
>>> import numpy as np
>>> from skcomplex.core.elements import Element, FieldElement
>>> from skcomplex.core.interactions import DiffusionInteraction
>>> from skcomplex.core.spaces import LatticeSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class CellElement(Element):
... def __init__(self, label, u: Annotated[float, StateSpec()] = 0.0):
... super().__init__(label)
... self.u = u
>>> space = LatticeSpace((40,), spacing=0.1, boundary='periodic')
>>> field = FieldElement('rod', element=CellElement, space=space, values={'u': np.eye(40)[20]})
>>> heat = DiffusionInteraction('heat', target='rod', state='u', coefficient=0.1)
>>> trajectory = Simulation(System(elements=[field], interactions=[heat]), dt=0.01).check().prepare()
>>> round(float(trajectory.run(until=1.0).trajectory_.values('rod.u')[-1].sum()), 6)
1.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'u',
coefficient: Callable[[SimpleNamespace], object] | float = 1.0,
) -> None:
super().__init__(label, target=target, state=state)
self.coefficient = coefficient
def _check(self) -> None:
"""Raise an error for a state that is not a name, or a coefficient that is neither callable nor a number.
Raises:
ValidationError:
The type declares states, the state is not a name, or the coefficient is neither callable nor a number.
"""
super()._check()
if not (callable(self.coefficient) or isinstance(self.coefficient, int | float)):
error_msg = f'The coefficient of {self.label!r} is {self.coefficient!r}. Give a function or a number.'
raise ValidationError(error_msg)
def _get_transport_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate `D laplacian(u)` of the state, the divergence of the diffusive flux.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespace of the field.
Returns:
laws:
The rate law of the state.
"""
coefficient = self.coefficient(namespace) if callable(self.coefficient) else self.coefficient
state = getattr(self._get_member(namespace, 'target'), self.state)
return [RateLaw(state, coefficient * laplacian(state))]
class RateInteraction(Interaction):
"""Rate `rate(ns)` and noise `noise(ns) dW` added to the highest time derivative of a state of a target.
On a state of order 1 the rate is its first derivative, and on a state of order `k` it adds to the `k`-th
derivative, so rates on one state superpose. With a divisor `d` the rate and the noise are divided by `d`, the name
of a state or parameter of the target, such as an inertia, a friction, or a capacitance, which backends can read.
A rate, a noise, or a divisor function takes the namespace of the interaction, with the namespaces `target` and
`sources` of the bound elements, the time `t`, and the symbols of the parameters a subclass declares. A number
stands for a constant, and a number on a vector state applies to every component. The noise is a Wiener process per
component and member by default, a process of a `WienerProcessSpec` to share or correlate it across states and
members, or a fractional Brownian motion of a Hurst index.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state the rate changes.
state:
The name of the state of the target that the rate changes.
rate:
A function of the namespace that returns the rate, or a number.
sources:
The paths of the elements or groups the rate and the noise read but do not change.
noise:
A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.
divisor:
The name of the state or parameter of the target that divides the rate and the noise, a function of the
namespace of the target that returns the divisor, such as a moment of inertia `m l^2`, or `None`.
process:
The Wiener processes that drive the noise, or `None` for one independent process per component and member.
index:
The process of `process` that drives the noise, or `None` for every process, one per component of a vector
state.
convention:
The stochastic calculus of the noise, `'ito'`, `'stratonovich'`, or `'anti-ito'`.
hurst:
The Hurst index in `(0, 1)` of a fractional Brownian motion that drives the noise, a number or a function
of the namespace that returns a parameter, or `None` for a Wiener process.
conserved:
Whether the noise is the divergence of a random flux on the links of the lattice of a field, so it keeps the
total of the state.
Raises:
ValidationError:
The state is not a name, a function is neither callable nor a number, the convention is unknown, the Hurst
index is outside `(0, 1)`, or the noise options cannot combine.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import RateInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> decay = RateInteraction('decay', target='stock', state='n', rate=lambda s: -0.5 * s.target.n)
>>> system = System(elements=[CompartmentElement('stock', n=8.0)], interactions=[decay])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=2.0).trajectory_
>>> round(float(trajectory.values('stock.n')[-1]), 3)
2.943
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'x',
rate: Callable[[SimpleNamespace], object] | float = 0.0,
sources: Annotated[Sequence[str], RoleSpec(variadic=True)] = (),
noise: Callable[[SimpleNamespace], object] | float | None = None,
divisor: str | Callable[[SimpleNamespace], object] | None = None,
process: Annotated[WienerProcessSpec | None, RoleSpec(noise=True)] = None,
index: int | None = None,
convention: str = 'ito',
hurst: Callable[[SimpleNamespace], object] | float | None = None,
conserved: bool = False,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.rate = rate
self.sources = sources
self.noise = noise
self.divisor = divisor
self.process = process
self.index = index
self.convention = convention
self.hurst = hurst
self.conserved = conserved
def _check_functions(self) -> None:
"""Raise an error for a state that is not a name, or a function that is neither callable nor a number."""
if not isinstance(self.state, str) or not self.state.isidentifier():
error_msg = (
f'{type(self).__name__} {self.label!r} changes the state {self.state!r}, which is not a name of a '
'state.'
)
raise ValidationError(error_msg)
for name in ('rate', 'noise'):
value = getattr(self, name)
if not (value is None or callable(value) or isinstance(value, int | float)):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace or a number.'
raise ValidationError(error_msg)
if not (
self.divisor is None
or callable(self.divisor)
or (isinstance(self.divisor, str) and self.divisor.isidentifier())
):
error_msg = (
f'The divisor of {self.label!r} is {self.divisor!r}. Give the name of a state or parameter, or a '
'function of the namespace of the target.'
)
raise ValidationError(error_msg)
def _check_noise(self) -> None:
"""Raise an error for an unknown convention, a Hurst index outside `(0, 1)`, or noise options that conflict."""
if self.convention not in CONVENTIONS:
error_msg = f'Unknown noise convention {self.convention!r}. Use one of {", ".join(CONVENTIONS)}.'
raise ValidationError(error_msg)
if not (
self.hurst is None or callable(self.hurst) or (isinstance(self.hurst, int | float) and 0 < self.hurst < 1)
):
error_msg = f'The Hurst index of {self.label!r} is {self.hurst!r}. Give a number in (0, 1) or a function.'
raise ValidationError(error_msg)
if self.index is not None and self.process is None:
error_msg = f'{self.label!r} picks the process {self.index!r} but binds no processes. Give `process`.'
raise ValidationError(error_msg)
options = [
name
for name, given in (('hurst', self.hurst), ('process', self.process), ('conserved', self.conserved))
if given
]
if len(options) > 1:
error_msg = f'The noise of {self.label!r} takes {" and ".join(options)}. Give at most one of them.'
raise ValidationError(error_msg)
def _check(self) -> None:
"""Raise an error for a state, functions, or noise options the interaction cannot take.
Raises:
ValidationError:
The type declares states, the state is not a name, a function is neither callable nor a number, the
convention is unknown, the Hurst index is outside `(0, 1)`, or the noise options cannot combine.
"""
super()._check()
self._check_functions()
self._check_noise()
def _evaluate(self, function: Callable[[SimpleNamespace], object] | float, namespace: SimpleNamespace) -> object:
"""Return the value of a function of the namespace, or a number as it is."""
return function(namespace) if callable(function) else function
def _broadcast(self, value: object, state: object) -> object:
"""Return a value with one component per component of a vector state, repeating a number."""
if isinstance(state, FreeModuleElement) and not isinstance(value, FreeModuleElement | list | tuple):
return vector(SR, [value] * len(state))
if isinstance(value, list | tuple):
return vector(SR, value)
return value
def _divide(self, value: object, namespace: SimpleNamespace) -> object:
"""Return a value divided by the divisor of the target, or the value without a divisor."""
if self.divisor is None:
return value
divisor = (
self.divisor(self._get_member(namespace, 'target'))
if callable(self.divisor)
else getattr(self._get_member(namespace, 'target'), self.divisor)
)
return operator.truediv(value, divisor)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the rate before the divisor, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.rate, namespace)
def _get_noise(self, namespace: SimpleNamespace) -> object:
"""Return the noise amplitude before the divisor, the hook a subclass overrides, or `None` for no noise."""
return None if self.noise is None else self._evaluate(self.noise, namespace)
def _get_noise_law(self, state: object, amplitude: object, namespace: SimpleNamespace) -> BasePart:
"""Return the noise law of an amplitude, by the process, the convention, and the Hurst index."""
if self.hurst is not None:
hurst = self._evaluate(self.hurst, namespace)
return FractionalNoiseLaw(state, amplitude, hurst=hurst, convention=self.convention)
process = None
if self.process is not None:
handle = self._get_member(namespace, 'process').dW
process = handle if self.index is None else handle[self.index]
return NoiseLaw(state, amplitude, convention=self.convention, process=process, conserved=self.conserved)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate law of the state and its noise law, both divided by the divisor.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the target and the sources.
Returns:
laws:
The rate law, and the noise law when the interaction has noise.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
rate = self._divide(self._broadcast(self._get_rate(namespace), state), namespace)
laws: list[BasePart] = [RateLaw(state, rate)]
noise = self._get_noise(namespace)
if noise is not None:
laws.append(self._get_noise_law(state, self._divide(noise, namespace), namespace))
return laws
Equations
gel.c[i]: diff(c_gel(t), t) == D_gel*laplacian(c_gel(t)) - partial(c_gel(t)*inverse_helmholtz(-partial(E_gel*partial(u_gel(t), 0) + T0_gel*c_gel(t)/(c_sat_gel + c_gel(t)), 0)/eta_gel, sqrt(gamma_gel/eta_gel)), 0)
gel.u[i]: diff(u_gel(t), t) == inverse_helmholtz(-partial(E_gel*partial(u_gel(t), 0) + T0_gel*c_gel(t)/(c_sat_gel + c_gel(t)), 0)/eta_gel, sqrt(gamma_gel/eta_gel))