BIO-LGCA

Last updated

In computational and mathematical biology, a biological lattice-gas cellular automaton (BIO-LGCA) is a discrete model for moving and interacting biological agents, [1] a type of cellular automaton. The BIO-LGCA is based on the lattice-gas cellular automaton (LGCA) model used in fluid dynamics. A BIO-LGCA model describes cells and other motile biological agents as point particles moving on a discrete lattice, thereby interacting with nearby particles. Contrary to classic cellular automaton models, particles in BIO-LGCA are defined by their position and velocity. This allows to model and analyze active fluids and collective migration mediated primarily through changes in momentum, rather than density. BIO-LGCA applications include cancer invasion [2] and cancer progression. [3]

Contents

Model definition

As are all cellular automaton models, a BIO-LGCA model is defined by a lattice , a state space , a neighborhood , and a rule . [4]

State space

The substructure of a BIO-LGCA lattice site with six velocity channels (corresponding to a 2D hexagonal lattice) and a single rest channel. In this case
b
=
6
{\displaystyle b=6}
,
a
=
1
{\displaystyle a=1}
, and the carrying capacity
K
=
7
{\displaystyle K=7}
. Channels 2, 3, 6 and 7 are occupied, thus the lattice configuration is
s
=
(
0
,
1
,
1
,
0
,
0
,
1
,
1
)
{\displaystyle \mathbf {s} =(0,1,1,0,0,1,1)}
, and the number of particles is
n
(
s
)
=
[?]
i
=
1
K
s
i
=
4
{\displaystyle n\left(\mathbf {s} \right)=\sum _{i=1}^{K}s_{i}=4}
. Hexnode.png
The substructure of a BIO-LGCA lattice site with six velocity channels (corresponding to a 2D hexagonal lattice) and a single rest channel. In this case , , and the carrying capacity . Channels 2, 3, 6 and 7 are occupied, thus the lattice configuration is , and the number of particles is .

For modeling particle velocities explicitly, lattice sites are assumed to have a specific substructure. Each lattice site is connected to its neighboring lattice sites through vectors called "velocity channels", , , where the number of velocity channels is equal to the number of nearest neighbors, and thus depends on the lattice geometry ( for a one-dimensional lattice, for a two-dimensional hexagonal lattice, and so on). In two dimensions, velocity channels are defined as . Additionally, an arbitrary number of so-called "rest channels" may be defined, such that , . A channel is said to be occupied if there is a particle in the lattice site with a velocity equal to the velocity channel. The occupation of channel is indicated by the occupation number . Typically, particles are assumed to obey an exclusion principle, such that no more than one particle may occupy a single velocity channel at a lattice site simultaneously. In this case, occupation numbers are Boolean variables, i.e. , and thus, every site has a maximum carrying capacity . Since the collection of all channel occupation numbers defines the number of particles and their velocities in each lattice site, the vector describes the state of a lattice site, and the state space is given by .

Rule and model dynamics

The states of every site in the lattice are updated synchronously in discrete time steps to simulate the model dynamics. The rule is divided into two steps. The probabilistic interaction step simulates particle interaction, while the deterministic transport step simulates particle movement.

Interaction step

Depending on the specific application, the interaction step may be composed of reaction and/or reorientation operators.

The reaction operator replaces the state of a node with a new state following a transition probability , which depends on the state of the neighboring lattice sites to simulate the influence of neighboring particles on the reactive process. The reaction operator does not conserve particle number, thus allowing to simulate birth and death of individuals. The reaction operator's transition probability is usually defined ad hoc form phenomenological observations.

The reorientation operator also replaces a state with a new state with probability . However, this operator conserves particle number and therefore only models changes in particle velocity by redistributing particles among velocity channels. The transition probability for this operator can be determined from statistical observations (by using the maximum caliber principle) or from known single-particle dynamics (using the discretized, steady-state angular probability distribution given by the Fokker-Planck equation associated to a Langevin equation describing the reorientation dynamics), [5] [6] and typically takes the form where is a normalization constant (also known as the partition function), is an energy-like function which particles will likely minimize when changing their direction of motion, is a free parameter inversely proportional to the randomness of particle reorientation (analogous to the inverse temperature in thermodynamics), and is a Kronecker delta which ensures that particle number before and after reorientation is unchanged.

The state resulting form applying the reaction and reorientation operator is known as the post-interaction configuration and denoted by .

Dynamics of the BIO-LGCA model. Every time step, the occupation numbers are changed stochastically by the reaction and/or reorientation operators in all lattice sites simultaneously during the interaction step. Subsequently, particles are deterministically moved to the same velocity channel on a neighboring node in the direction of their velocity channel, during the transport step. Colors in the sketch are used to track the dynamics of the particles of individual nodes. This sketch assumes a particle-conserving rule (no reaction operator). Hexdynamics.png
Dynamics of the BIO-LGCA model. Every time step, the occupation numbers are changed stochastically by the reaction and/or reorientation operators in all lattice sites simultaneously during the interaction step. Subsequently, particles are deterministically moved to the same velocity channel on a neighboring node in the direction of their velocity channel, during the transport step. Colors in the sketch are used to track the dynamics of the particles of individual nodes. This sketch assumes a particle-conserving rule (no reaction operator).

Transport step

After the interaction step, the deterministic transport step is applied synchronously to all lattice sites. The transport step simulates the movement of agents according to their velocity, due to the self-propulsion of living organisms.

During this step, the occupation numbers of post-interaction states will be defined as the new occupation states of the same channel of the neighboring lattice site in the direction of the velocity channel, i.e. .

A new time step begins when both interaction and transport steps have occurred. Therefore, the dynamics of the BIO-LGCA can be summarized as the stochastic finite-difference microdynamical equation

Example interaction dynamics

A hexagonal BIO-LGCA model of polar swarming. In this model, cells preferentially change their velocities to be parallel to the neighborhood's momentum. Lattice sites are colored according to their orientation, following the color wheel. Empty sites are white. Periodic boundary conditions were used.

The transition probability for the reaction and/or reorientation operator must be defined to appropriately simulate the modeled system. Some elementary interactions and the corresponding transition probabilities are listed below.

Random walk

In the absence of any external or internal stimuli, cells may move randomly without any directional preference. In this case, the reorientation operator may be defined through a transition probability

A hexagonal BIO-LGCA model of excitable media. In this model, the reaction operator favors the rapid reproduction of particles within velocity channels, and the slow death of particles within rest channels. Particles in rest channels inhibit the reproduction of particles in velocity channels. The reorientation operator is the random walk operator in the text. Lattice sites are brightly colored the more motile particles are present. Resting particles are not shown. Periodic boundary conditions were used.

where . Such transition probability allows any post-reorientation configuration with the same number of particles as the pre-reorientation configuration , to be picked uniformly.

Simple birth and death process

If organisms reproduce and die independently of other individuals (with the exception of the finite carrying capacity), then a simple birth/death process can be simulated [3] with a transition probability given by where , are constant birth and death probabilities, respectively, is the Kronecker delta which ensures only one birth/death event happens every time step, and is the Heaviside function, which makes sure particle numbers are positive and bounded by the carrying capacity .

A square BIO-LGCA model of cells interacting adhesively. Cells move preferentially in the direction of the cell density gradient. Lattice sites are colored with increasingly darker blue colors with increasing cell density. Empty nodes are colored white.Periodic boundary conditions are used.

Adhesive interactions

Cells may adhere to one another by cadherin molecules on the cell surface. Cadherin interactions allow cells to form aggregates. The formation of cell aggregates via adhesive biomolecules can be modeled [7] by a reorientation operator with transition probabilities defined as

A square BIO-LGCA model of cells indirectly interacting chemotactically. In this model, cells produce a diffusing chemoattractant with a certain half-life. Cells preferentially move in the direction of the chemoattractant gradient. Lattice sites are additively colored with a darker blue tint with increasing cell density, and with a darker yellow tint with increasing chemoattractant concentration. Empty lattice sites are colored white. Periodic boundary conditions were used.

where is a vector pointing in the direction of maximum cell density, defined as , where is the configuration of the lattice site within the neighborhood , and is the momentum of the post-reorientation configuration, defined as . This transition probability favors post-reorientation configurations with cells moving towards the cell density gradient.

Mathematical analysis

Since an exact treatment of a stochastic agent-based model quickly becomes unfeasible due to high-order correlations between all agents, [8] the general method of analyzing a BIO-LGCA model is to cast it into an approximate, deterministic finite difference equation (FDE) describing the mean dynamics of the population, then performing the mathematical analysis of this approximate model, and comparing the results to the original BIO-LGCA model.

First, the expected value of the microdynamical equation is obtained where denotes the expected value, and is the expected value of the -th channel occupation number of the lattice site at at time step . However, the term on the right, is highly nonlinear on the occupation numbers of both the lattice site and the lattice sites within the interaction neighborhood , due to the form of the transition probability and the statistics of particle placement within velocity channels (for example, arising from an exclusion principle imposed on channel occupations). This non-linearity would result in high-order correlations and moments among all channel occupations involved. Instead, a mean-field approximation is usually assumed, wherein all correlations and high order moments are neglected, such that direct particle-particle interactions are substituted by interactions with the respective expected values. In other words, if are random variables, and is a function, then under this approximation. Thus, we can simplify the equation to where is a nonlinear function of the expected lattice site configuration and the expected neighborhood configuration dependent on the transition probabilities and in-node particle statistics.

From this nonlinear FDE, one may identify several homogeneous steady states, or constants independent of and which are solutions to the FDE. To study the stability conditions of these steady states and the pattern formation potential of the model, a linear stability analysis can be performed. To do so, the nonlinear FDE is linearized as where denotes the homogeneous steady state , and a von Neumann neighborhood was assumed. In order to cast it into a more familiar finite difference equation with temporal increments only, a discrete Fourier transform can be applied on both sides of the equation. After applying the shift theorem and isolating the term with a temporal increment on the left, one obtains the lattice-Boltzmann equation [4] where is the imaginary unit, is the size of the lattice along one dimension, is the Fourier wave number, and denotes the discrete Fourier transform. In matrix notation, this equation is simplified to , where the matrix is called the Boltzmann propagator and is defined as The eigenvalues of the Boltzmann propagator dictate the stability properties of the steady state: [4]

Applications

Constructing a BIO-LGCA for the study of biological phenomena mainly involves defining appropriate transition probabilities for the interaction operator, though precise definitions of the state space (to consider several cellular phenotypes, for example), boundary conditions (for modeling phenomena in confined conditions), neighborhood (to match experimental interaction ranges quantitatively), and carrying capacity (to simulate crowding effects for given cell sizes) may be important for specific applications. While the distribution of the reorientation operator can be obtained through the aforementioned statistical and biophysical methods, the distribution of the reaction operators can be estimated from the statistics of in vitro experiments, for example. [9]

BIO-LGCA models have been used to study several cellular, biophysical and medical phenomena. Some examples include:

Related Research Articles

Continuum mechanics is a branch of mechanics that deals with the deformation of and transmission of forces through materials modeled as a continuous medium rather than as discrete particles. The French mathematician Augustin-Louis Cauchy was the first to formulate such models in the 19th century.

A phonon is a collective excitation in a periodic, elastic arrangement of atoms or molecules in condensed matter, specifically in solids and some liquids. A type of quasiparticle in physics, a phonon is an excited state in the quantum mechanical quantization of the modes of vibrations for elastic structures of interacting particles. Phonons can be thought of as quantized sound waves, similar to photons as quantized light waves.

<span class="mw-page-title-main">Fokker–Planck equation</span> Partial differential equation

In statistical mechanics and information theory, the Fokker–Planck equation is a partial differential equation that describes the time evolution of the probability density function of the velocity of a particle under the influence of drag forces and random forces, as in Brownian motion. The equation can be generalized to other observables as well. The Fokker-Planck equation has multiple applications in information theory, graph theory, data science, finance, economics etc.

<span class="mw-page-title-main">Hamiltonian mechanics</span> Formulation of classical mechanics using momenta

In physics, Hamiltonian mechanics is a reformulation of Lagrangian mechanics that emerged in 1833. Introduced by Sir William Rowan Hamilton, Hamiltonian mechanics replaces (generalized) velocities used in Lagrangian mechanics with (generalized) momenta. Both theories provide interpretations of classical mechanics and describe the same physical phenomena.

<span class="mw-page-title-main">Mie scattering</span> Scattering of an electromagnetic plane wave by a sphere

In electromagnetism, the Mie solution to Maxwell's equations describes the scattering of an electromagnetic plane wave by a homogeneous sphere. The solution takes the form of an infinite series of spherical multipole partial waves. It is named after German physicist Gustav Mie.

<span class="mw-page-title-main">Path integral formulation</span> Formulation of quantum mechanics

The path integral formulation is a description in quantum mechanics that generalizes the stationary action principle of classical mechanics. It replaces the classical notion of a single, unique classical trajectory for a system with a sum, or functional integral, over an infinity of quantum-mechanically possible trajectories to compute a quantum amplitude.

In physics, the S-matrix or scattering matrix relates the initial state and the final state of a physical system undergoing a scattering process. It is used in quantum mechanics, scattering theory and quantum field theory (QFT).

In physics, a parity transformation is the flip in the sign of one spatial coordinate. In three dimensions, it can also refer to the simultaneous flip in the sign of all three spatial coordinates :

<span class="mw-page-title-main">LSZ reduction formula</span> Connection between correlation functions and the S-matrix

In quantum field theory, the Lehmann–Symanzik–Zimmermann (LSZ) reduction formula is a method to calculate S-matrix elements from the time-ordered correlation functions of a quantum field theory. It is a step of the path that starts from the Lagrangian of some quantum field theory and leads to prediction of measurable quantities. It is named after the three German physicists Harry Lehmann, Kurt Symanzik and Wolfhart Zimmermann.

<span class="mw-page-title-main">Mathematical formulation of the Standard Model</span> Mathematics of a particle physics model

This article describes the mathematics of the Standard Model of particle physics, a gauge quantum field theory containing the internal symmetries of the unitary product group SU(3) × SU(2) × U(1). The theory is commonly viewed as describing the fundamental set of particles – the leptons, quarks, gauge bosons and the Higgs boson.

In quantum mechanics, the probability current is a mathematical quantity describing the flow of probability. Specifically, if one thinks of probability as a heterogeneous fluid, then the probability current is the rate of flow of this fluid. It is a real vector that changes with space and time. Probability currents are analogous to mass currents in hydrodynamics and electric currents in electromagnetism. As in those fields, the probability current is related to the probability density function via a continuity equation. The probability current is invariant under gauge transformation.

Ewald summation, named after Paul Peter Ewald, is a method for computing long-range interactions in periodic systems. It was first developed as the method for calculating the electrostatic energies of ionic crystals, and is now commonly used for calculating long-range interactions in computational chemistry. Ewald summation is a special case of the Poisson summation formula, replacing the summation of interaction energies in real space with an equivalent summation in Fourier space. In this method, the long-range interaction is divided into two parts: a short-range contribution, and a long-range contribution which does not have a singularity. The short-range contribution is calculated in real space, whereas the long-range contribution is calculated using a Fourier transform. The advantage of this method is the rapid convergence of the energy compared with that of a direct summation. This means that the method has high accuracy and reasonable speed when computing long-range interactions, and it is thus the de facto standard method for calculating long-range interactions in periodic systems. The method requires charge neutrality of the molecular system to accurately calculate the total Coulombic interaction. A study of the truncation errors introduced in the energy and force calculations of disordered point-charge systems is provided by Kolafa and Perram.

In condensed matter physics and crystallography, the static structure factor is a mathematical description of how a material scatters incident radiation. The structure factor is a critical tool in the interpretation of scattering patterns obtained in X-ray, electron and neutron diffraction experiments.

<span class="mw-page-title-main">Random phase approximation</span>

The random phase approximation (RPA) is an approximation method in condensed matter physics and in nuclear physics. It was first introduced by David Bohm and David Pines as an important result in a series of seminal papers of 1952 and 1953. For decades physicists had been trying to incorporate the effect of microscopic quantum mechanical interactions between electrons in the theory of matter. Bohm and Pines' RPA accounts for the weak screened Coulomb interaction and is commonly used for describing the dynamic linear electronic response of electron systems. It was further developed to the relativistic form (RRPA) by solving the Dirac equation.

Static force fields are fields, such as a simple electric, magnetic or gravitational fields, that exist without excitations. The most common approximation method that physicists use for scattering calculations can be interpreted as static forces arising from the interactions between two bodies mediated by virtual particles, particles that exist for only a short time determined by the uncertainty principle. The virtual particles, also known as force carriers, are bosons, with different bosons associated with each force.

Heat transfer physics describes the kinetics of energy storage, transport, and energy transformation by principal energy carriers: phonons, electrons, fluid particles, and photons. Heat is thermal energy stored in temperature-dependent motion of particles including electrons, atomic nuclei, individual atoms, and molecules. Heat is transferred to and from matter by the principal energy carriers. The state of energy stored within matter, or transported by the carriers, is described by a combination of classical and quantum statistical mechanics. The energy is different made (converted) among various carriers. The heat transfer processes are governed by the rates at which various related physical phenomena occur, such as the rate of particle collisions in classical mechanics. These various states and kinetics determine the heat transfer, i.e., the net rate of energy storage or transport. Governing these process from the atomic level to macroscale are the laws of thermodynamics, including conservation of energy.

The semiconductor luminescence equations (SLEs) describe luminescence of semiconductors resulting from spontaneous recombination of electronic excitations, producing a flux of spontaneously emitted light. This description established the first step toward semiconductor quantum optics because the SLEs simultaneously includes the quantized light–matter interaction and the Coulomb-interaction coupling among electronic excitations within a semiconductor. The SLEs are one of the most accurate methods to describe light emission in semiconductors and they are suited for a systematic modeling of semiconductor emission ranging from excitonic luminescence to lasers.

Lagrangian field theory is a formalism in classical field theory. It is the field-theoretic analogue of Lagrangian mechanics. Lagrangian mechanics is used to analyze the motion of a system of discrete particles each with a finite number of degrees of freedom. Lagrangian field theory applies to continua and fields, which have an infinite number of degrees of freedom.

In machine learning, the kernel embedding of distributions comprises a class of nonparametric methods in which a probability distribution is represented as an element of a reproducing kernel Hilbert space (RKHS). A generalization of the individual data-point feature mapping done in classical kernel methods, the embedding of distributions into infinite-dimensional feature spaces can preserve all of the statistical features of arbitrary distributions, while allowing one to compare and manipulate distributions using Hilbert space operations such as inner products, distances, projections, linear transformations, and spectral analysis. This learning framework is very general and can be applied to distributions over any space on which a sensible kernel function may be defined. For example, various kernels have been proposed for learning from data which are: vectors in , discrete classes/categories, strings, graphs/networks, images, time series, manifolds, dynamical systems, and other structured objects. The theory behind kernel embeddings of distributions has been primarily developed by Alex Smola, Le Song , Arthur Gretton, and Bernhard Schölkopf. A review of recent works on kernel embedding of distributions can be found in.

The Peierls substitution method, named after the original work by Rudolf Peierls is a widely employed approximation for describing tightly-bound electrons in the presence of a slowly varying magnetic vector potential.

References

  1. Deutsch, Andreas; Nava-Sedeño, Josué Manik; Syga, Simon; Hatzikirou, Haralampos (2021-06-15). "BIO-LGCA: A cellular automaton modelling class for analysing collective cell migration". PLOS Computational Biology. 17 (6): e1009066. Bibcode:2021PLSCB..17E9066D. doi: 10.1371/journal.pcbi.1009066 . ISSN   1553-7358. PMC   8232544 . PMID   34129639.
  2. Reher, David; Klink, Barbara; Deutsch, Andreas; Voss-Böhme, Anja (2017-08-11). "Cell adhesion heterogeneity reinforces tumour cell dissemination: novel insights from a mathematical model". Biology Direct. 12 (1): 18. doi: 10.1186/s13062-017-0188-z . ISSN   1745-6150. PMC   5553611 . PMID   28800767.
  3. 1 2 Böttger, Katrin; Hatzikirou, Haralambos; Voss-Böhme, Anja; Cavalcanti-Adam, Elisabetta Ada; Herrero, Miguel A.; Deutsch, Andreas (2015-09-03). Alber, Mark S (ed.). "An Emerging Allee Effect Is Critical for Tumor Initiation and Persistence". PLOS Computational Biology. 11 (9): e1004366. Bibcode:2015PLSCB..11E4366B. doi: 10.1371/journal.pcbi.1004366 . ISSN   1553-7358. PMC   4559422 . PMID   26335202.
  4. 1 2 3 "Mathematical Modeling of Biological Pattern Formation", Cellular Automaton Modeling of Biological Pattern Formation, Modeling and Simulation in Science, Engineering and Technology, Boston, MA: Birkhäuser Boston, pp. 45–56, 2005, doi:10.1007/0-8176-4415-6_3, ISBN   978-0-8176-4281-5 , retrieved 2021-05-25
  5. Nava-Sedeño, J. M.; Hatzikirou, H.; Peruani, F.; Deutsch, A. (2017-02-27). "Extracting cellular automaton rules from physical Langevin equation models for single and collective cell migration". Journal of Mathematical Biology. 75 (5): 1075–1100. doi:10.1007/s00285-017-1106-9. ISSN   0303-6812. PMID   28243720. S2CID   32456636.
  6. Nava-Sedeño, J. M.; Hatzikirou, H.; Klages, R.; Deutsch, A. (2017-12-05). "Cellular automaton models for time-correlated random walks: derivation and analysis". Scientific Reports. 7 (1): 16952. arXiv: 1802.04201 . Bibcode:2017NatSR...716952N. doi:10.1038/s41598-017-17317-x. ISSN   2045-2322. PMC   5717221 . PMID   29209065.
  7. Bussemaker, Harmen J. (1996-02-01). "Analysis of a pattern-forming lattice-gas automaton: Mean-field theory and beyond". Physical Review E. 53 (2): 1644–1661. Bibcode:1996PhRvE..53.1644B. doi:10.1103/physreve.53.1644. ISSN   1063-651X. PMID   9964425.
  8. Ovaskainen, Otso; Somervuo, Panu; Finkelshtein, Dmitri (2020-10-28). "A general mathematical method for predicting spatio-temporal correlations emerging from agent-based models". Journal of the Royal Society Interface. 17 (171): 20200655. doi:10.1098/rsif.2020.0655. PMC   7653394 . PMID   33109018.
  9. Dirkse, Anne; Golebiewska, Anna; Buder, Thomas; Nazarov, Petr V.; Muller, Arnaud; Poovathingal, Suresh; Brons, Nicolaas H. C.; Leite, Sonia; Sauvageot, Nicolas; Sarkisjan, Dzjemma; Seyfrid, Mathieu (2019-04-16). "Stem cell-associated heterogeneity in Glioblastoma results from intrinsic tumor plasticity shaped by the microenvironment". Nature Communications. 10 (1): 1787. Bibcode:2019NatCo..10.1787D. doi:10.1038/s41467-019-09853-z. ISSN   2041-1723. PMC   6467886 . PMID   30992437.
  10. Mente, Carsten; Prade, Ina; Brusch, Lutz; Breier, Georg; Deutsch, Andreas (2010-10-01). "Parameter estimation with a novel gradient-based optimization method for biological lattice-gas cellular automaton models". Journal of Mathematical Biology. 63 (1): 173–200. doi:10.1007/s00285-010-0366-4. ISSN   0303-6812. PMID   20886214. S2CID   12404555.
  11. Bussemaker, Harmen J.; Deutsch, Andreas; Geigant, Edith (1997-06-30). "Mean-Field Analysis of a Dynamical Phase Transition in a Cellular Automaton Model for Collective Motion". Physical Review Letters. 78 (26): 5018–5021. arXiv: physics/9706008 . Bibcode:1997PhRvL..78.5018B. doi:10.1103/PhysRevLett.78.5018. ISSN   0031-9007. S2CID   45979152.
  12. Fuks, Henryk; Lawniczak, Anna T. (2001). "Individual-based lattice model for spatial spread of epidemics". Discrete Dynamics in Nature and Society. 6 (3): 191–200. doi: 10.1155/s1026022601000206 . hdl: 1807/82157 .
  13. Ilina, Olga; Gritsenko, Pavlo G.; Syga, Simon; Lippoldt, Jürgen; La Porta, Caterina A. M.; Chepizhko, Oleksandr; Grosser, Steffen; Vullings, Manon; Bakker, Gert-Jan; Starruß, Jörn; Bult, Peter (2020-08-24). "Cell–cell adhesion and 3D matrix confinement determine jamming transitions in breast cancer invasion". Nature Cell Biology. 22 (9): 1103–1115. doi:10.1038/s41556-020-0552-6. ISSN   1476-4679. PMC   7502685 . PMID   32839548.